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I.  INTRODUCTION 


A.  BACKGROUND 

The  finite  element  (FEM)  method  is  the  representation  of  structures  by  a 
discretized  set  of  elements  and  nodes  for  the  purpose  of  engineering 
calculations.  This  method  has  permeated  almost  every  form  of  engineering  due 
to  its  adaptability  and  capability  to  handle  everything  from  heat  transfer  to 
structural  dynamics.  The  FEM  utilizes  shape  functions  and  boundary  conditions 
formulated  from  the  prediction  of  how  an  element  reacts  to  the  stimuli  provided 
and  interpolated  throughout  each  element.  As  this  thesis  focuses  on  the  field  of 
structural  dynamics,  the  finite  element  method  is  typically  paired  with  the  theory 
of  elasticity  in  order  to  formulate  the  shape  functions  required  for  calculation. 

The  theory  of  elasticity  is  a  mathematical  interpretation  of  the  deformation 
of  structures  and  objects  [1].  For  the  purposes  of  structural  dynamics  and  its 
applicability  to  this  thesis,  we  will  use  the  linear  formulation,  which  assumes  the 
following:  The  structure  experiences  small  deformations  (infinitesimal  strain),  the 
material  has  a  linear  stress/strain  relationship,  and  the  material  does  not 
experience  yielding  (as  yielding  is  a  non-linear  process).  The  theory  of  elasticity 
can  be  applied  using  any  coordinate  system  or  in  the  direct  tensor  form.  The 
theory  of  elasticity  is  the  basis  of  many  common  engineering  equations  such  as 
Flooke’s  law  and  the  bi-harmonic  wave  equation;  however,  we  will  focus  on  a 
simplified  version  of  the  theory  called  the  Euler-Bernoulli  beam  theory. 

The  Euler-Bernoulli  beam  theory,  as  will  be  re-iterated  in  Chapter  II,  was 
derived  from  the  theory  of  elasticity  based  on  the  load  carrying  and  deflection 
characteristics  of  beams  [1 , 2].  To  simplify  the  calculation  of  stress  and  deflection 
of  a  beam,  the  Euler-Bernoulli  beam  theory  assumes  the  neutral  axis  does  not 
change  length,  all  cross  sections  remain  plane  to  the  neutral  axis,  and  the 
material  is  homogeneous  between  ends  as  demonstrated  in  Figure  1 . 
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Cross-section 


Figure  1.  Euler-Bernouili  Beam  Assumption  Characteristics.  Source:  [3]. 

As  described  earlier,  the  finite  element  method  and  Euler-Bernouili  beam 
theory  are  combined  in  structural  dynamics  to  predict  the  response  of  a  beam 
under  static  and  dynamic  forces  [4],  The  Euler-Bernouili  beam  theory  and  finite 
element  method  are  applicable  to  a  three-dimensional  structure;  however,  we  will 
simplify  their  application  to  two  dimensions  for  this  thesis  based  on  the 
assumption  that  movement  will  be  planar.  The  simulated  model  of  the  beam 
comprises  several  discrete  beam  elements  assembled  computationally  in 
accordance  with  the  finite  element  method,  whereas  each  element  of  the  beam 
reflects  and  responds  as  dictated  by  the  Euler-Bernouili  beam  theory.  While  there 
is  no  guidance  as  to  how  many  elements  need  to  exist,  the  larger  number  of 
elements  that  comprise  the  model  will  develop  a  better  prediction  for  the 
response  of  a  prototype  beam.  However,  as  the  number  of  elements  increase,  so 
does  the  amount  of  calculations.  Other  considerations  into  how  many  elements 
should  comprise  the  simulation  may  be  influenced  by  testing  equipment  and 
measurement  capabilities.  As  such,  the  resolution  of  the  elements  is  just  one 
possible  cause  of  which  the  simulated  model  may  be  required  to  be  updated  to 
better  reflect  a  prototype  beam. 

Of  the  many  issues  that  exist  between  simulation  and  prototype,  one  of 
the  most  common  is  an  inaccurate  beam  dynamic  response  prediction  as  shown 
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in  [2,  5].  There  are  many  potential  causes  for  the  error,  including  material 
property  inaccuracies,  manufacturing  discrepancies,  measurement  errors, 
improper  assumptions,  etc.  Conversely,  some  of  these  errors  may  also  indicate 
locations  of  potential  damage  in  the  prototype,  if  the  simulated  model  proved  to 
be  an  accurate  representation  upon  initial  build.  One  of  the  indicators  of  a 
discrepancy  between  simulation  and  prototype  is  the  beam  harmonic  response  in 
response  to  a  force  input,  also  known  as  the  frequency  response.  The  frequency 
response  function  (FRF)  is  a  measure  of  the  magnitude  and  phase  of  beam 
displacement  per  unit  force  as  described  in  [6].  An  example  of  this  can  be  found 
in  Figure  2. 


Image  created  using  [7], 

Figure  2.  Frequency  Response  Function 
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The  FRF  provides  details  as  to  resonance  frequencies  for  the  structure, 
which  correlate  to  the  modes  of  bending  and  correlate  to  the  associated 
resonance  frequencies  calculated  from  the  finite  element  model  (FEM).  The  FRF 
can  be  computationally  derived  or  physically  measured.  The  computation  is 
described  more  thoroughly  in  Chapter  II,  but  consists  of  determining  the  beam 
response  at  a  specific  node  and  specific  degree  of  freedom  given  a  harmonic 
forcing  function  applied  at  another  degree  of  freedom  in  the  beam.  This 
computation  can  also  be  found  in  [4],  Similarly,  the  measured  FRF  is  the  ratio  of 
the  harmonic  excitation  force  provided  by  an  impact  hammer  or  shaker,  to  the 
response  force  measured  by  an  accelerometer,  as  demonstrated  by  [7],  Many 
types  of  programs  exist  that  can  gather  and  calculate  FRF  data  from  such  an 
experiment,  but  typically  are  limited  by  equipment  costs  and  in  frequency 
bandwidth.  By  comparing  the  FRF  from  the  computation  to  the  measured  FRF, 
we  are  able  to  determine  if  there  are  discrepancies  between  the  two.  If  there  are, 
we  can  use  a  method  known  as  “Sensitivity-Based  Finite  Element  Model 
Updating  and  Structural  Damage  Detection”  [5,  8]  to  determine  where  the 
discrepancies  for  the  FRF’s  originate 

The  sensitivity  matrix  used  for  updating  the  model  is  calculated  from  the 
FEM  eigenvectors  and  the  FEM  reaction  to  a  perturbation  at  one  element  [5]. 
The  process  of  this  is  discussed  further  in  Chapter  II,  but  the  resulting  matrix 
relates  the  change  in  resonance  frequencies  given  a  change  in  a  physical 
parameter  of  an  element.  However,  this  method  is  restricted  to  the  amount  of 
modes/resonance  frequencies  that  are  shown  in  both  measured  and  calculated 
FRF  plots,  and  thus  the  bandwidth  of  the  experimental  equipment.  Typically, 
equipment  falls  far  short  of  capturing  all  modes  of  vibration  [5,  9],  which  can  lead 
to  mode  truncation  errors  when  comparing  the  FRFs  and  give  incomplete  or 
incorrect  indications  of  errors.  It  is  from  this  issue  that  this  thesis  will  focus  and 
develop  a  new  method  for  model  updating  and  damage  detection. 

The  method  that  this  thesis  explores  involves  the  installation  of  either  a 
massless  spring  or  a  lumped  mass  into  the  structure  in  an  attempt  to  provide  new 
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and  unique  data  for  damage  detection  or  model  updating.  Frequency  domain 
structural  synthesis  as  shown  in  [5,  10]  is  an  exact  method  to  calculate  changes 
to  the  FRF  of  a  structure  after  structural  modification  or  substructure  coupling. 
The  synthesis  technique  is  useful  because  of  its  capability  to  handle  large  model 
reductions,  frequency  dependencies,  and  provide  an  exact  solution  based  on  the 
baseline  response  of  the  structure  and  the  impedance  descriptions  of  the 
structural  modifications.  Previous  explorations  in  this  field,  such  as  [11],  utilized 
the  artificial  synthesis  of  pin  connections  to  provide  unique  data  and  proved 
successful. 

B.  LITERATURE  REVIEW 

“Classic  beam  theory,”  as  derived  from  the  theory  of  elasticity,  was 
developed  by  Jacob  Bernoulli  and  applied  by  Daniel  Bernoulli  and  Leonhard 
Euler  to  describe  static  and  dynamic  beam  elements  in  response  from  stimuli. 
This  is  more  thoroughly  described  in  [12].  The  dynamic  response  produces  a 
FRF  from  which  the  resonance  frequencies  can  be  determined  and  later  used  for 
damage  detection  or  model  updating.  As  discussed  earlier,  the  dynamic 
response  can  be  determined  two  ways,  though  measurements  or  though 
calculation.  The  measurement  process  consists  of  placing  one  or  more 
accelerometers  along  the  prototype  beam  and  exciting  the  structure  with  either  a 
special  hammer,  or  attaching  an  electrodynamic  shaker.  Software  is  used 
provide  a  suitable  excitation  signal  for  an  electrodynamic  shaker  and  to  measure 
the  dynamic  response  of  the  structure  to  produce  a  highly  detailed  FRF,  as 
demonstrated  in  [7],  However,  the  measureable  bandwidth  captures  far  fewer 
modes  than  the  full  dynamic  response  of  the  beam,  estimated  on  the  order  of 
1 01— 1 02  times  less,  as  described  in  [5].  The  calculated  dynamic  response  is 
determined  by  calculating  the  eigenvectors  and  eigenvalues  of  the  global  mass 
and  stiffness  matrices  developed  from  the  FEM  [5].  The  eigenvalues  directly 
correlate  to  the  resonance  frequencies  in  the  structure,  and  the  eigenvalues  are 
the  mass-normalized  positions  of  the  nodes  of  the  FEM  for  each  resonance 
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frequency,  otherwise  known  as  a  mode.  Further  explanation  can  be  found  in  [13, 
4],  This  method  can  determine  as  many  frequencies  and  modes  as  there  are 
degrees  of  freedom  (DOF)  in  the  model,  but  higher  frequencies  will  have 
increasingly  larger  error  from  the  actual  frequencies  unless  node  count  increases 
to  infinite,  which  is  obviously  not  possible.  As  such,  one  possible  measure  for 
determining  a  minimum  amount  of  nodes  is  the  convergence  of  the  resonance 
frequencies  within  the  measurable  bandwidth. 

Another  beam  theory  exists  known  as  the  Timoshenko  theory  and  is 
explored  in  [1,  12],  As  Timoshenko  and  Flan  explain,  instead  of  assuming  that 
cross-sections  remain  perpendicular  to  the  neutral  axis,  the  Timoshenko  theory 
assumes  there  is  angular  deformation,  but  the  neutral  axis  remains  the  same 
length.  This  method  has  proven  to  be  more  accurate  in  the  representation  of  the 
bending  of  structures  where  the  cross-section  does  not  remain  perpendicular  to 
the  neutral  axis,  however  has  minor  gains  in  accuracy  for  a  rigid,  smaller 
structure  [12],  As  such,  this  thesis  will  continue  to  use  the  Euler-Bernoulli 
beam  theory. 

As  discussed  earlier,  the  current  effort  is  to  utilize  both  the  measured  FRF 
and  the  calculated  FRF  to  determine  the  possible  location  of  either  damage  or 
necessity  for  model  updating,  which  are  described  in  [5,  14],  When  calculating 
the  sensitivity  matrix,  the  number  of  modes  that  are  within  the  measureable 
bandwidth  typically  are  fewer  than  the  number  of  modes  available  to  the  whole 
FEM.  This  creates  a  rectangular  matrix  that  is  typically  not  full  rank,  and  thus 
computationally  cannot  be  reliably  solved  for  damage  detection  or  for  model 
updating.  In  an  effort  to  correct  this,  several  methods  have  been  attempted. 

One  method  to  provide  new  information  into  the  measured  FRF  is  to 

modify  the  physical  structure  with  known  boundary  conditions,  as  demonstrated 

by  [15].  This  method  would  physically  place  a  pin  condition  or  add  a  mass  at  a 

node  of  the  prototype  beam  to  impart  a  “perturbed  boundary  condition”  (PBC). 

Each  modification  would  alter  the  FRF  and  provide  more,  relevant  data  to  the 

calculations.  Flowever,  the  implementation  of  the  method  has  in  increasing  risk  of 
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introducing  inaccuracies  in  the  characterization  of  boundary  conditions. 
Additionally,  this  method  requires  increasing  amounts  of  testing  and  a  laboratory 
work,  becoming  increasingly  expensive. 

Another  proposed  method  for  model  updating  uses  the  anti-resonance 
frequencies  of  the  FRF  as  opposed  to  the  resonance  frequencies,  as  explored  in 
[16].  Where  resonance  in  an  FRF  is  a  peak  value  where  a  node  has  the  most 
displacement  at  a  specific  mode,  anti-resonance  is  where  there  is  minimal  to  no 
movement  at  a  specific  node  given  a  certain  frequency  on  the  FRF.  Typical 
accelerometers  send  a  milli-voltage  back  to  the  measuring  equipment  based  on 
the  rate  of  movement  they  experience,  demonstrated  by  [7],  Because  the  anti¬ 
resonance  point  is  determined  to  have  no  movement,  the  signal  the 
accelerometers  provide  would  be  subject  to  noise,  making  the  anti-resonance 
point  difficult  to  measure  accurately. 

Finally,  this  thesis  focuses  on  modifying  a  method  that  introduces  “artificial 
boundary  conditions”  to  provide  new  and  relevant  information  to  the  sensitivity 
matrix  as  originally  explored  by  [5].  The  method  computationally  modifies  test 
data  through  a  process  called  synthesis  to  apply  a  boundary  condition,  thus  no- 
longer  requiring  the  expensive  and  time  consuming  physical  alterations. 
Boundary  conditions  can  be  systematically  applied  to  any  FEM  and  quickly: 
calculate  the  eigenvectors  and  eigenvalues,  calculate  a  new  FRF,  create  a  new 
sensitivity  matrix  based  on  the  new  conditions,  and  provide  relevant  and  new 
information  for  model  updating  and  damage  detection.  The  synthesized  test  data 
compared  to  the  modified  FEM  will  provide  new  information  to  the  modal 
sensitivity  matrix  of  which,  with  the  previous  calculations  appended  to  the  new 
information,  provides  a  much  more  reasonable  and  reliable  method  for  model 
updating  and  damage  detection. 

Based  on  the  ABC  method  developed  in  [5,  11],  an  updated  procedure 

was  developed  to  identify  parameter-specific  optimal  artificial  boundary 

conditions  [8].  When  calculating  the  differences  in  frequencies  of  a  damaged  or 

undamaged  beam,  the  likelihood  to  correctly  find  discrepancies  between  the 
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model  and  the  prototype  per  each  FEM  element  is  predominantly  determined  by 
the  decomposition  of  the  sensitivity  matrix.  In  this  method,  artificial  boundary 
conditions  are  chosen  based  on  the  QR  method  for  matrix  decomposition.  The 
matrix  decomposition  is  accomplished  by  using  MATLAB  proprietary  coding  to 
provide  a  stable  decomposition,  explored  in  [8,  17].  The  method  then  analyzes 
the  sequence  of  pivots  of  the  decomposition  and  utilizes  the  rows  most  likely  to 
find  damage  within  each  element.  The  method  eventually  creates  a  composite 
square,  full  rank  sensitivity  matrix  that  is  far  more  effective  and  accurate  in 
determining  damage  or  updating  the  model. 

C.  SCOPE  OF  THESIS 

This  thesis  will  explore  the  use  of  an  artificial  spring  or  dimensionless 
mass  added  to  the  beam  and  connected  to  ground  in-lieu  of  a  pin  connection.  By 
varying  the  spring  stiffness  at  a  single  node,  we  determine  if  new  and  relevant 
frequency  data  and  subsequent  sensitivity  matrices  can  be  used  in  damage 
detection/model  updating.  Additionally,  the  beam  should  react  the  same  as  a  pin 
connection  when  using  a  spring  of  near-infinite  stiffness. 

For  the  purpose  of  the  exploration  of  this  new  method,  finite  element 
model  will  be  limited  to  a  two  dimensional  Euler-Bernoulli  beam  structure  and 
discretized  using  the  finite  element  method.  This  model  reflects  a  prototype 
structure  for  practical  comparison.  If  there  is  new  and  relevant  information,  this 
new  method  explores  the  accuracy  of  its  application. 

D.  NOTES  ON  THESIS 

As  mentioned  earlier,  the  FEM  reflects  a  prototype  beam.  Beam 
measurements  are  conducted  using  Pulse  Reflex  software  [7]  and  calculations 
are  conducted  in  MATLAB  software  [17]. 
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II.  THEORY 


A.  STRUCTURAL  DYNAMICS 

It  is  theorized  a  structure  can  experience  a  near  infinite  number  of 
resonance  frequencies  as  a  continuous  beam  has  an  infinite  number  of  degrees 
of  freedom(DOF).  However,  while  we  can  use  a  significant  number  of  elements, 
the  finite  element  process  discretizes  the  beam  and  number  of  degrees  of 
freedom  to  form  a  finite  element  model  based  on  the  Euler-Bernoulli  beam 
theory. 

The  Euler-Bernoulli  Beam  describes  the  relationship  of  a  beam  in  reaction 
to  static  and  dynamic  forces  derived  through  Newton’s  second  law  and  Rayleigh 
energy  equations  as  found  in  [4],  As  described  before,  this  beam  theory  operates 
on  the  following  assumptions: 

1 .  The  cross-section  plane  of  the  beam  remains  plane  during 
deformation. 

2.  The  beam  has  a  neutral  axis  that  does  not  experience  expansion  or 
contraction. 

3.  Cross  sections  that  are  perpendicular  to  the  neutral  axis  remain 
perpendicular  through  deformation. 

4.  The  beam  is  homogenous  and  has  a  linear  modulus  of  elasticity 
E(x). 

5.  Primary  stresses  perpendicular  to  the  neutral  axis  are  negligible 
when  compared  to  primary  stresses  parallel  with  the  neutral  axis. 

6.  The  rotary  inertia  of  the  beam  is  neglected  in  the  moment  equation. 

7.  Displacements  are  relatively  small. 

Figure  3  displays  the  coordinate  system  of  a  typical  2-D  Euler-Bernoulli 

Beam. 
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With  these  assumptions,  the  differential  equation  of  motion  for  a 
transverse  vibration  of  a  beam  is  [4]: 


8 2 
dx2 


f 

El 

V 


dV 

dx 2  y 


+  pA 


f(x,t)  — >  0  <  x  <  L  , 


(2.1) 


where  El  is  the  flexural  rigidity  consisting  of  the  modulus  of  elasticity  (E)  and  the 
moment  of  inertia(l)  for  the  cross  sectional  area,  p  is  the  mass  density,  A  is  the 
cross  sectional  area,  f  is  the  force  experiences  by  the  beam.  However,  when  the 
beam  is  experiencing  transverse  loads  only,  and  no  forcing  function,  the 
equilibrium  equation,  as  explained  by  [4],  reduces  to: 


d2 

dx2 
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El 
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dV 

dx2  y 


fix,  t) 


(2.2) 


1.  Shape  Functions 

We  approximate  the  dynamic  transverse  displacement  of  the  beam  over 
the  length  of  the  beam  using  the  following  equation,  as  explained  by  [4],: 


4 

y(x,t)  =  Yji//i(x)yi(t)  (2.3) 

i=i 

Where  ip  represents  the  shape  functions  (assumed  modes)  of  the  beam. 
Because  the  general  solution  to  equation  2.2  is: 


10 


(2.4) 


y(x)  =  c1+c2 


fx3 

+  Co 

2 

+  C4 

^X1 

J 

4 

as  explained  by  [4], we  can  assume  that  the  shape  functions  take  the  form  of  a 
cubic.  However,  those  shape  functions  must  fulfill  the  boundary  conditions  found 
in  equation  (2.5)  and  displayed  in  Figure  4  [4]: 


1. ^(0)  =  1  y/  '(0)  =  y/(L)  =  y/  '(L)  =  0 

2.  y/  '(0)  =  1  ^(0)  =  ^(L)  =  yr  '(L)  =  0 

3.  ^(L)  =  1  ^(0)  =  y/  '(0)  =  y/  '(L)  =  0 

4.  y/  '(L)  =  1  ^(0)  =  y/  '(0)  =  y/(L)  =  0 


Figure  4.  Cubic  Hermite  Shape  Functions.  Source:  [18]. 

By  entering  the  boundary  conditions  from  equation  (2.5)  into  (2.4),  we  get 
the  following  shape  function  for  the  beam  as  developed  by  [4]: 
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(2.6) 


By  using  the  equations  for  strain  and  kinetic  energy  in  a  beam,  we  are 
able  to  determine  the  coefficients  for  the  stiffness  or  mass  matrices  using  the 
following  equations  as  identified  in  [4]: 

L  M  |? 

k,j  =  | EIy/i  y/j  dx  i,j=l,2,3,4  (2.7) 

0 

L 

mtj  =  | pAy/i  y j  dx  i,j=l,2,3,4  (2.8) 

0 

By  plugging  in  the  shape  functions  (2.6)  into  the  coefficient  equations 
(2. 7, 2. 8),  we  are  able  to  determine  the  stiffness  and  mass  matrices  of  a  uniform 
beam  element  that  is  loaded  only  by  end  shears  and  moments.  The  result  is  the 
following,  as  calculated  in  [4]: 
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(2.9) 


(2.10) 
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Equations  2.9  and  2.10  are  used  in  the  construction  of  the  finite  element 


model. 


2.  Finite  Element  Model 

The  FEM  is  a  structure  that  is  comprised  of  multiple  Euler-Bernoulli  beam 
elements  connected  in  an  effort  to  predict  overall  dynamic  response.  As  such, 
the  governing  equation  for  a  dynamic  system  is  the  following,  as  derived  by  [4]: 


MM+MW+MM ={/<*>} 


(2.11) 


In  this  instance,^  and  ^  are  the  mass  and  stiffness  matrices  for  the 

[cl 

whole  system,  otherwise  known  as  global  matrices.1  J  is  the  global  damping 
[M]  j  [K] ,  ancj  [C]  each  js  an  nxn  matrix  where  n  is  the  number 


coefficient  matrix. 


of  degrees  of  freedom  in  the  structure.  is  an  nxl  array  comprising  of  the 


l'Xj  is  the 


external  forces  on  the  structure  at  all  degrees  of  freedom. 

Ui 

displacement  and  rotation  position  of  the  beam,  1  >  is  the  transverse  and 

rotational  velocity,  and  1  ’  is  the  transverse  and  rotational  acceleration  of  the 
beam. 


For  the  purposes  of  this  thesis,  we  will  assume  that  our  beam  experiences 
undamped,  free  vibrations.  The  reason  we  assume  the  structure  is  undamped  is 
because  damping  cannot  be  predicted,  and  has  little  influence  on  the  frequency 
of  resonance  peaks.  Typically,  modal  damping  is  measured  in  a  vibration  test; 
however,  without  a  test,  damping  is  assumed  at  2%  (f;=0.02)  and  proportional  to 
the  stiffness  and  mass  matrices.  The  influence  of  damping  is  determined  from 
the  typical  conversion  from  an  undamped  resonance  frequency  to  a  damped 
resonance  frequency,  found  in  [4]: 


co. 


=  co„ 
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(2.12) 


Based  on  equation  (2.12),  the  typical  assumed  damping  would  change  the 
resonance  frequency  by  0.04%,  which  is  a  small  enough  change  to  safely 
neglect. 

By  neglecting  damping  in  the  system,  the  governing  equation  for  a 
dynamic  system  changes  to  the  following,  as  derived  by  [4]: 


(2.13) 


The  assembly  of  global  [m]  and  [k]\s  determined  by  the  degrees  of 

freedom  and  node  connections.  Because  we  are  using  a  four  DOF  beam 
element,  the  assembly  of  these  matrices  follows  this  pattern: 
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(2.14) 


(2.15) 


The  final  product  produces  a  set  of  global  matrices  that  are  nxn  where 
n=number  of  DOFs=(DOF/node)*m  elements. 


B.  NATURAL  FREQUENCY  SENSITIVITY 


One  of  the  benefits  of  the  FEM  is  the  capability  to  determine  the  modal 

frequencies  and  mode  shapes  inherent  to  the  structure  without  the  necessity  of 

performing  a  test.  Flowever,  after  building  and  testing  a  prototype,  we  want  to  be 

able  to  use  the  differences  between  the  FEM  and  the  prototype  modal 

frequencies  to  find  errors  in  the  FEM.  Some  possibilities  of  errors  could  be 
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incorrect  assumptions  in  the  model  or  unmodeled  joint  mechanics  in  a  prototype. 
A  method  of  determining  these  errors  is  the  use  of  the  Sensitivity  matrix 

1.  Modal  Frequencies 

Continuing  with  the  assumption  of  an  undamped  bar  experiencing  free 
vibrations,  we  can  assume  a  harmonic  solution  for  equation  (2.13)  is  the 
following,  as  determined  by  [4,  9]: 


M  =  {(/>\CejM ,  (2.16) 

where  <j>  is  the  eigenvector  for  the  mode  of  vibration.  The  eigenvector  is  a  product 
of  the  eigenvalue  problem  that  details  the  position  of  each  degree  of  freedom  in 
the  structure  for  the  corresponding  eigenvalue.  The  resulting  algebraic 

eigenvalue  problem  as  described  in  [4],  [14],  and  [19]  is: 

[[*MM]fo}=0  (2.17) 

We  are  looking  for  more  than  a  trivial  solution  of  {4>}  ^  0  and  therefore  the 

eigenvalue,^  ,  must  take  on  a  value  to  render  the  matrix  non-invertable  and  thus 

yielding  non  trivial  solutions  for {^} .  As  such,  the  values  of  \  are  found  by  setting 
the  determinant  of  the  difference  to  zero,  as  determined  by  [4,  9]: 

Det\[K]-\  [M]]  =  0  (2.18) 

There  are  as  may  roots  to  this  equation  as  there  are  DOF  in  the  structure. 
As  \  =o)f,  the  eigenvalue  can  quickly  be  converted  to  the  resonance 

frequencies  of  the  structure.  This  method  is  validated  by  the  use  of  a  secondary 
dynamic  system  defined  though  analytical  and  omitted  coordinates. 
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2.  Analytical  and  Omitted  Coordinates 

As  alluded  in  Chapter  I,  it  is  theorized  that  any  given  structure  may  have 
resonance  frequencies  on  the  order  of  1 02— 1 03  ,  which  is  further  discussed  in  [5]. 
However,  a  vibration  test  may  only  measure  modal  parameters  on  the  order  of 
101.  In  order  to  validate  the  applicability  of  the  FEM  and  the  vibration  test,  we 
take  a  look  at  the  analytical  and  omitted  coordinates  of  the  beam.  The  analytical 
coordinates  (ASET)  are  the  measured  coordinates  where  instrumentation  is 
applied  to  the  structure  and  is  annotated  in  equations  using  “a.”  The  omitted 
coordinates  (OSET)  are  the  coordinates  that  are  not  directly  measured  and  are 
annotated  in  equations  using  “o.”  For  the  purposes  of  the  FEM,  the  OSET  is  of 
finite  dimension;  however,  the  OSET  is  of  infinite  dimension  in  a  vibration  test. 

The  OSET  is  also  referred  as  a  secondary  dynamic  system  and  is  defined 
by  the  process  of  representing  a  large  order  dynamic  system  with  a  smaller  order 
dynamic  system  [14].  The  secondary  dynamic  system  is  obtained  by  fully 
constraining  a  coordinate  from  the  ASET.  This  effectively  removes  a  row  and 
column  from  the  stiffness  and  mass  matrices  in  the  FEM.  This  introduces 
singularities  into  the  undamped  system  which  translate  to  modal  frequencies. 

Re-evaluating  equation  (2.17)  for  a  forced  system  of  analytical  and 
omitted  coordinates,  as  described  in  [14,  20]: 
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(2.19) 


However,  we  assume  there  is  no  excitation  acting  on  the  omitted 
coordinates.  As  such,  equation  (2.19)  becomes,  [14,  20]: 
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(2.20) 


If  we  examine  the  two  row  partitions  of  equation  2.20,  we  get  the  following, 
[14,  20]: 
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+  +  *„}]  =  /.  (2-21) 

[^]{^}  +  [««,]{^}-®2[[^]{^}  +  [M„]{^}]=0  (2.22) 

Using  equation  (2.22),  we  can  find  the  exact  dynamic  relationship  between 
the  omitted  and  analytic  coordinates  as  calculated  in  [14]: 

W  =  [I -fflXK]'1  [-K^Koa  +G>2K£Moa'\{xa}  (2.23) 

While  equation  (2.23)  is  an  exact  reduction  of  the  model,  it  provides  some 
issues  with  regards  to  the  conditioning  of  the  matrix.  As  such,  we  use  the 
Improved  Reduction  System  (IRS)  developed  and  tested  in  [8,  14,  20].  Using  the 
static  Test  Analysis  Model  (TAM)  reduction: 

R}=®2RHM»f  K„,]R)  (2.24) 

The  new  exact  model  reduction  is  frequency  independent  [8, 14,  20]: 

{*>}  =  [-K0oK„„  +  (-/r>oa  -  KooMUoKoo Koa ) M-sl,Kstat ] (xa }  (2.25) 

However,  in  order  for  equation  (2.23)  to  work,  [i  J'  needs  to  be 

invertible.  A  matrix  can  be  inverted  by  dividing  the  adjoint  (Adj)  by  the 
determinant  (Det).  In  equation  form,  the  formula  for  an  inversion  is  [8,  14,  20].: 


w1 


Det[A] 


Adj  [A] 


(2.26) 


So,  using  equation  (2.23),  the  following  needs  to  be  satisfied  [8,  14,  20]. 


J‘  =  r  1  M[l-r»X'Mm]  (2.27) 

L  J  Det\\  -ruAjM,,,  J  L  J 

Similarly  to  what  was  identified  in  equation  (2.18),  the  roots  of  the 
following  characteristic  equation  [8,  14,  20].: 


Det[l-^XxXo]  '=0 
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(2.28) 


are  the  resonance  frequencies  for  the  OSET  system,  made  up  of  [£oo]and[Moo], 
when  the  ASET  frequencies  are  constrained  to  ground. 

3.  Building  the  Modal  Sensitivity  Matrix 

The  purpose  of  the  modal  sensitivity  matrix  is  to  assist  in  predicting  the 
change  in  dynamic  response  of  a  structure  based  on  the  change  of  a  given 
parameter.  As  such,  it  is  developed  by  perturbing,  or  increasing,  a  parameter  a 
minor  amount,  such  as  1  percent,  on  a  single  element  of  the  structure  and 
determining  change  in  the  resonance  frequency  for  that  mode.  To  put  it 
mathematically,  as  found  in  [14,  21]: 


S 


ij 


rU, 

dPj 


(2.29) 


Where  \  is  the  eigenvalue  for  mode  i,  Pj  is  the  perturbed  parameter  at 

element  j,  and  Stj  is  the  sensitivity  value  for  mode  /  and  with  a  modification  in 
element  j. 

The  sensitivity  matrix  is  calculated  by  perturbing  an  element  in  the 
structure  and  calculating  the  modal  response.  More  specifically,  we  calculate  the 
change  in  the  resonance  frequencies  by  using  the  eigenvector  for  that  mode,  as 
demonstrated  in  [4,  21]. 


8K  .  SM 

'  —  A; 


dp,  '  dpf 


Ml  ’ 


(2.30) 


where  {$}  is  the  eigenvector  for  mode  /,  {$}r  is  the  transpose  of  the  eigenvector, 


and 


cK  dM 
dPj  '  dPj 


is  the  change  in  the  stiffness  and  mass  matrices  with  respect  to 


the  change  in  a  parameter  at  element  j.  {4}  is  of  size  nxl  where  n  is  the  number 


18 


of  modes/DOFs  in  the  FEM,  and 


is  sized  n  x  n.  Equation  2.30  yields 


8K  t  8M 

dpj  '  8Pj 

a  singular  entry  of  the  sensitivity  matrix 

The  sensitivity  matrix  should  be  sized  n  x  m,  where  m  is  the  number  of 
elements  in  the  structure,  and  n  is  the  number  of  modes.  For  example: 
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The  importance  of  the  size  of  the  sensitivity  matrix  will  be  discussed  later 
in  this  chapter.  Several  different  types  of  parameters  can  vary  in  the  structure. 
Some  examples  of  these  are  the  flexural  rigidity  (E*l),  width,  or  depth  of  the 
beam.  Once  the  sensitivity  value  has  been  calculated,  it  can  be  used  to  predict 
future  eigenvalues/resonance  frequencies  given  an  arbitrary  change  in  that 
parameter.  As  found  in  [14,  21]: 


A*  =4  +Sij-(p*-pj),  (2.32) 

where  A*  is  the  new  eigenvalue/resonance  frequency  given  the  updated 
parameter  p*  at  element  j.  Alternatively,  we  can  reduce  Equation  2.19  to  simply 

needing  to  know  the  change  in  a  resonance  frequency  based  on  the  change  of  a 
parameter  at  a  given  element.  [21] 


=  Stj  •  (APj )  (2.33) 

For  example,  assuming  we  need  to  know  what  the  change  is  to  the  4th 
mode  of  vibration  due  to  a  change  in  the  El  value  of  element  2,  the  calculation 
would  look  something  like  this: 


A  a>l  =  AT4  -S42-  (A  El) 
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(2.34) 


4.  Determining  the  Modal  Sensitivity  Matrix  Using  Eigenvectors 
and  Eigen  Values 

The  development  of  equation  (2.29)  as  developed  by  [14]  and  [22]  starts 
with  equation  (2.17)  where  [m]  and  [.fifjare  the  mass  and  stiffness  matrices  as 

developed  by  the  FEM.  As  identified  in  equation  (2.30)  the  parameter  that  is 
being  modified  may  be  a  part  of  either  or  both  matrices.  For  example,  if  modal 

parameter  Pj  is  the  density  of  the  jth  element,  only  the  mass  will  be  effected. 

However,  if  Pj  is  the  thickness  of  the  jth  element,  Pj  is  then  a  function  of  the 

stiffness  and  mass  matrices.  As  such,  the  derivative  for  both  matrices  must  be 
calculated,  as  developed  by  [9,  13,  21, 23]: 


8K 

dP.i 


(2.34) 


dM 
dPj _ 


(2.35) 


As  such,  if  we  take  the  derivative  of  equation  (2.17)  with  respect  to  pjt  we 
get  the  following  equation  [9,  13,  21 , 23]: 


dK 

dPj 


dPj  dpj 


=  0 


(2.36) 


By  pre-multiplying  equation  (2.36)  by  {$}r,  we  get  the  following  as 
demonstrated  by  [14]  and  [22]: 


5  dM  d\ 

A - -M 

dpj  dp 


=  0 


(2.37) 


By  recognizing  [m]  =  [m]7  and  [7T]  =  [AT]7  and  invoking  equation  (2.17) 


gives  the  following: 
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T 


(2.38) 


6K  „  0M  61  ,, 

- X - -M 

dp  j  dp.  dp . 


=  0 


If  we  assume  a  mass  normalized  mode  shape{$}r[M]{$}  =  1,  we  get  the 
following: 


04 

dp, 


dK 

dp, 


A 


dM 

dp, 


(2.39) 


Which  is  equation  (2.29).  However,  if  we  are  conducting  a  test  on  a 
prototype,  we  most  likely  will  not  know  the  changed  parameter,  or  on  which 
element  has  been  altered.  As  such,  we  can  conduct  a  vibration  test  of  the 
structure  and  determine  the  resonance  frequencies  from  the  frequency  response 
functions. 


C.  FREQUENCY  RESPONSE  FUNCTION 

The  frequency  response  function  (FRF)  is  a  measure  of  the  dynamic 
response  of  the  structure  and  can  be  both  calculated  and  measured.  The  FRF  is 
a  transfer  function  from  the  input  (force)  to  the  output  (acceleration).  The 
derivation  of  the  calculation  can  be  found  in  [4], 

The  quantity  [[A'j-QfM]]  as  discussed  when  determining  the  modal 

frequencies  is  known  as  the  impedance,  [Z] .  The  impedance  of  a  structure  is  a 
calculation  of  the  structure’s  resistance  to  dynamic  movement  when  subjected  to 
harmonic  motion.  The  inverse  of  the  impedance  is  the  FRF,  [h].  With  regard  to 

the  viewpoint  of  a  transfer  function,  the  input  is  a  harmonic  signal,  at  some 
frequency  Q,  and  the  output  is  the  beam  response  [21]: 

[[ATl-ntAf]]-1  =[z]-!  =[^T]  (2.40) 
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The  calculated  FRF,  [h],  records  the  beam  response  of  every  DOF  as 

each  node  is  separately  excited,  for  a  specific  excitation  frequency.  As  such,  if 
the  beam  response  was  needed  from  1  hertz  to  1000  hertz  with  a  0.1  hertz  step, 
the  FRF  calculation,  equation  (2.40),  would  need  to  be  performed  for  each  step 
of  the  range. 

Once  the  FRF  matrix  has  been  calculated,  we  can  extract  the 
displacement  of  the  DOF  by  multiplying  the  FRF  by  the  force  at  that  frequency, 
as  explained  in  [21]: 


{X(0)}=[ff]{F(n)}  (2.41) 

By  expanding  this  to  an  n  DOF  system,  the  system  of  equations  would 
look  something  like  this: 
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(2.42) 


Where  Xn  is  the  complex  amplitude  response  of  the  n-th  DOF,  Hjk  is  the  j- 

th  DOF  response  due  to  harmonic  excitation  at  the  k-th  DOF,  and  Fn  if  the 
harmonic  force  at  the  n-th  DOF.  All  of  these  terms  are  frequency  dependent. 

When  measuring  the  beam  response  for  the  FRF,  high  fidelity 
measurement  devices  are  required  to  determine  the  structure’s  dynamic 
response  to  excitation.  As  will  be  discussed  in  Chapter  III,  each  node  needs  to  be 
excited  separately  while  measuring  the  response  of  the  beam  with  measurement 
software  for  processing  to  develop  a  full  FRF  matrix.  Ultimately,  a  measured  FRF 
will  produce  a  set  of  resonance  frequencies,  go,  that  can  be  converted  into 
eigenvalues.  This  allows  for  the  most  complete  comparison  when  performing 
damage  detection/model  updating  calculations. 
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D.  DAMAGE  DETECTION  AND  MODEL  UPDATING 

When  we  are  looking  at  the  differences  between  a  FEM  and  a  structure 
using  the  FRF  comparison  method,  we  are  typically  looking  for  “hidden”  damage, 
or  something  that  would  produce  relatively  minor  differences  in  resonance 
frequencies  between  the  model  and  the  structure.  Model  updating  is  the  process 
of  updating  the  FEM  based  on  the  dynamic  response  of  a  recently  built 
prototype.  Damage  detection  is  the  process  of  identifying  problematic  elements 
based  on  the  differences  between  dynamic  response  of  the  updated  FEM  and 
the  response  of  the  same  structure  that  has  many  hours  of  use.  In  both  cases, 
the  goal  of  using  this  method  is  to  identify  where  the  FEM  does  not  correlate  to 
the  reaction  of  the  prototype. 

Ultimately,  we  want  to  use  the  resonance  frequency  data  for  damage 
detection  and/or  model  updating.  For  the  purposes  of  this  thesis,  we  will  assume 
the  parameter  in  question  is  the  flexural  rigidity  (E*l),  and  the  sensitivity  matrix 
has  been  calculated  to  reflect  as  such.  To  accomplish  damage  detection  or 
model  updating,  we  change  the  perspective  of  how  we  view  equation  (2.33)  to 
the  following  equations,  as  detailed  in  [14]: 


AA^SAEL 

(2.43) 

(\fEM  ~  Aj  Measured  )  =  ^ ijAEI j 

(2.44) 

{A/t}  =  [5]{A£/} 

(2.45) 

[/>]  1  {A/L}  =  {A£7} 

(2.46) 

Where  [S]  is  the  full  sensitivity  matrix,  {AA}  is  a  vector  consisting  of  the 

difference  in  eigenvalues  for  each  mode,  and  {A£7}  is  a  vector  identifying  which 

element  the  suspected  parameter  differs.  However,  in  order  for  this  process  to 
work,  the  sensitivity  matrix  must  be  invertible  and  full  rank.  This  means  that  the 
measuring  software  will  need  to  capture  as  many  modes  in  the  structure  as  there 
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are  elements,  at  a  minimum.  Unfortunately,  the  bandwidth  of  most  measuring 
equipment  is  not  capable  of  measuring  higher  order  frequencies  [5]  and  thus 
cannot  include  modes  above  some  practical  bandwidth  in  the  calculations.  This 
limitation  leads  to  a  non-square  sensitivity  matrix,  also  known  as  modal 
truncation.  The  other  limitation  of  the  sensitivity  matrix  being  full  rank  becomes 
an  issue  when  you  have  a  symmetrical  structure.  The  rank  refers  to  the  number 
of  linearly  independent  rows  and  columns  in  a  matrix.  If  the  structure  is 
symmetrical,  the  columns  of  the  sensitivity  matrix  may  be  close  enough  in  value 
to  be  considered  linearly  dependent.  In  an  effort  to  correct  these  issues,  we 
implement  a  derivative  of  the  Artificial  Boundary  Condition  (ABC)  method 
detailed  in  [5]. 

E.  STRUCTURAL  SYNTHESIS 

A  typical  structural  dynamic  design  cycle  will  involve  repeated  creation  of 
finite  element  models  due  to  unacceptable  errors  or  results  from  the 
previous  iteration  [21].  For  large  FEMs,  this  can  be  very  computationally 
demanding  and  expensive  in  manpower  and  time.  The  synthesis  method 
eliminates  the  necessity  to  re-create  and  re-solve  an  entire  FEM  after  the 
baseline  evaluation.  The  method  is  able  to  accomplish  this  by  applying  models 
of  structural  changes  to  the  baseline  dynamic  response  and  effectively 
transform  the  baseline  response  into  the  modified  response.  The  speed  of  the 
computation  comes  into  effect  by  only  including  the  FRF  for  the  DOF  involved 
in  the  synthesis  and  additional  DOF  in  which  we  have  interest. 

DOFs  involved  in  synthesis  typically  fall  under  two  conditions:  The 
connection  set  (CSET),  or  the  internal  set  (ISET).  The  CSET  describes  the  DOFs 
of  connection  of  a  modification  to  the  original  structure.  In  this  instance,  the 
CSET  will  be  limited  to  a  single  DOF  for  a  boundary  condition  or  a  spring/mass 
addition.  DOFs  that  might  not  be  involved  in  the  connection,  but  need  to  be 
recorded,  are  known  as  the  internal  set  (ISET).  Together,  these  points  are  known 
as  the  analytical  set  (ASET).  For  example,  if  we  have  a  beam  modification  that 
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will  be  connected  to  nodes  6  and  10,  and  we  want  to  record  the  translational 
response  at  node  1 5  and  1 6,  we  identify  DOF  11,12,19  and  20  as  the  CSET  DOF 
and  DOF  29  and  31  as  the  ISET  DOF.  This  example  is  displayed  in  Figure  5  and 
equation  (2.47) 


Beam  Modification 

Figure  5.  16-Element  Structure  with  Beam  Modification. 


CSET  =  [11,12,19, 20] 

ISET  =  [29,31]  (2.47) 

ASET  =  [CSET,  ISET]  =  [1 1,12,1 9, 20, 29, 3 1] 

Once  these  sets  have  been  identified,  we  can  pull  them  out  of  the  primary 
FRF  matrix: 
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Or,  with  the  analytical  set  definitions,  we  can  re-write  equation  (2.47)  to 
something  like  this,  as  demonstrated  in  [21]: 
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In  the  process  of  determining  the  effects  of  the  modification,  we  need  to 
determine  the  force  exerted  by  the  modification  on  the  structure.  As  such,  we 
modify  equation  (2.41)  to  solve  for  force  as  opposed  to  displacement: 
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{f;}=-[ak+®=am]{x;}=-[az]{x;}  (2.50) 

Where  {f,*}  is  the  force  of  the  modification  at  the  CSET  DOF,  jx*}  is  the 

response  the  modification  “sees,”  and  [AX  +  <y2AM]or[AZ]  is  the  added 
modification.  Additionally,  if  there  is  any  external  force  on  the  modification  or 
structure,  it  can  be  annotated  with  {fcw}.  Therefore,  as  described  in  [5],  the 
overall  force  applied  to  the  CSET  is: 

(Fc  }  =  {f«}  +  {f;}  =  {f“}-[AZ]{x;J  (2.51) 

As  such,  we  can  expand  equation  (2.49)  to  the  following: 
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(2.52) 


Based  on  our  definition  of  ASET,  we  can  re-write  equation  (2.52)  to  the 
following: 
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CC 

ac 


(2.53) 


We  can  break  equation  (2.53)  into  its  system  of  equations  like  so: 


{x;}=[ff.,]{F;}  +  [ffJ{F“}  (2.54) 

(2-55) 

By  using  equation  (2.45)  and  solving  for  jxj}  in  equation  (2.49),  we  get 
the  following: 


{ K  }  =  [/  +  [He e  ]  [ AZ  ( o)]}'  [Hc a  ]  { } 


(2.56) 
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Now,  by  using  equation  (2.50)  and  equation  (2.56)  in  equation  (2.55),  we 
get  the  following: 


The  final  form  of  the  synthesized  matrix  is: 


(2.57) 


{X.'HCKM  (2-58) 

As  such,  the  synthesized  FRF  matrix  is: 


(2.59) 


1.  Artificial  Boundary  Condition  Synthesis 

The  Artificial  Boundary  Condition  (ABC)  method  computationally  applies 
restrictions  to  the  test  data  to  simulate  a  pin  connected  to  ground  [5].  Finding 
these  resonance  frequencies  for  known  changes  helps  to  provide  new 
information  to  the  sensitivity  matrix  and  provide  reliability  in  damage  detection. 
For  the  purpose  of  this  thesis,  we  will  be  identifying  a  connection  point  to  apply 
our  conditions,  and  we  will  be  interested  in  the  response  of  the  remaining  DOFs 
of  the  structure.  As  such,  [CSET]=restricted  DOF  and  [ISET]=remaining  DOF.  An 
example  of  an  artificial  pin  constraint  is  displayed  in  Figure  6. 


Figure  6.  1 6-Element  Structure  with  Artificial  Boundary  Condition. 


Because  the  CSET  will  have  no  movement,  {Xc}  =  {0} .  As  such,  equation 
(2.48)  changes  to  the  following: 
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(2.60) 


By  using  the  top  row  of  equation  (2.60),  we  can  solve  for  the  forces  from 
synthesis  (the  asterisk  indicating  the  effects  of  the  modification): 


By  plugging  equation  (2.61)  into  the  bottom  row  of  equation  (2.60)  and  re¬ 
arranging  terms,  we  get  the  following: 

(2-62) 

Similarly  to  equation  (2.58),  jx*}  represents  the  beam  response  as  a 
result  of  the  ABC  synthesis.  Therefore,  the  modified  FRF  response  is: 


[«;]  =  [[«, (2-63) 

Similarly  to  when  we  analyzed  the  analytical  and  omitted  coordinates,  the 
act  of  putting  a  pin  connection  into  a  FEM  structure  removes  the  column  and  row 
associated  with  the  DOF  in  the  stiffness  and  mass  matrices.  This  an  easy 
method  to  check  the  synthesis  calculations.  By  recalculating  the  baseline  FRF 
matrix  with  the  associated  modifications,  we  can  overlay  the  graph  onto  the 
synthesized  version.  If  done  correctly,  the  check  and  the  synthesized  FRF  should 
be  exact. 


2.  Artificial  Spring  or  Mass  Synthesis 

The  focus  of  this  thesis  is  the  introduction  of  an  artificial,  massless  spring 
or  lumped  mass  into  the  FEM  and  determining  the  response.  An  example  of 
using  synthesizing  an  artificial  spring  is  in  Figure  7. 
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Figure  7.  1 6-Element  Structure  with  Artificial  Spring. 


The  two  DOF  massless  spring  is  connected  from  the  node  to  ground. 
Similar  to  the  ABC,  the  DOF  that  would  have  been  introduced  from  the  end 
connected  to  ground  is  no  longer  considered.  The  spring  value  is  synthesized 
into  the  structure  FRF  without  adding  DOF  to  the  system. 

For  the  instance  of  a  lumped  mass,  we  are  attaching  the  mass  directly  to 
the  node.  Because  it  is  an  idealized  condition  known  as  “lumped,”  there  are  no 
rotational  inertial  components  to  account  for  and  so  the  added  artificial  mass  will 
only  effect  the  translational  DOF. 

For  both  instances,  the  CSET  is  the  DOF  that  we  are  synthesizing  the 
spring  or  mass  onto,  and  the  ISET  are  the  remaining  DOFs  of  the  structure. 

Because  we  are  using  only  a  single  value  to  synthesize,  equation  (2.58) 
becomes: 


As  demonstrated  in  [21].  In  this  instance: 


(2.64) 


[AZ(®)]  =  A  kspring 

or  (2.65) 

[_AZ(co)]  =  -co2Ammass 

The  benefit  of  this  method  is  the  wide  range  of  variability  for  either  value. 
By  using  a  mass,  we  should  see  the  FRF  resonance  peaks  shift  down  in  the 
measurable  bandwidth.  By  using  a  spring,  we  should  see  the  FRF  resonance 
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values  shift  upward  in  the  measureable  bandwidth.  Ultimately,  a  spring  of  infinite 
stiffness  will  act  as  a  pin  connection. 

F.  COMPOSITE  SENSITIVITY  MATRIX 

As  discussed  earlier,  one  of  the  primary  issues  in  damage  detection  and 
model  updating  is  modal  truncation.  Modal  truncation  occurs  due  to  a  difference 
in  the  number  of  modes  that  can  be  identified  in  a  vibration  test  as  opposed  to 
those  that  can  be  calculated  from  a  FEM.  For  instance,  a  2-dimensional  beam 
with  20  elements  will  produce  40  modes;  however,  testing  equipment  may  be 
limited  to  measuring  up  to  1000  Hz,  and  thus  can  only  capture  14  modes.  As 
such,  we  calculate  the  sensitivity  matrix  with  the  14  modes  and  get  a  sensitivity 
matrix  that  is  14  modes  x  20  elements.  In  this  configuration,  the  sensitivity  matrix 
is  underdetermined  and  rank  deficient.  As  such,  it  cannot  reliably  be  used  in 
damage  detection  or  model  updating  as  its  use  in  equation  (2.46)  will  produce 
false  indications  of  damage  in  an  attempt  to  solve  the  equation. 

Another  issue  presents  itself  when  you  are  given  a  symmetric  structure 
such  as  a  free-free  beam.  If  20  modes  where  available  for  our  20  element  beam, 
we  have  a  fully  developed  sensitivity  matrix,  but  we  may  still  have  an  issue  as  the 
structure  is  symmetrical.  Symmetry  in  a  structure  will  produce  a  symmetric 
sensitivity  matrix,  thus  making  it  rank  deficient  since  it  will  not  have  enough 
unique  values  to  solve  the  system  of  equations. 

Because  of  modal  truncation  and/or  rank  deficiency  in  the  sensitivity 
matrix,  we  use  the  Artificial  Boundary  Condition  method  to  calculate  a  new 
sensitivity  matrix  and  new  frequencies.  The  process  of  which  requires  the  FEM  to 
be  re-calculated  for  each  condition  presented. 

Additionally,  as  this  is  for  the  initial  development  of  the  sensitivity  matrix, 
we  create  a  second  FEM  of  the  structure  identical  except  for  a  known  amount  of 
damage.  This  allows  us  to  develop  a  {A/l}  and  {A£7}  matrix  based  on  the 
conditions  provided. 


30 


For  the  development  of  the  composite  sensitivity  matrix,  we  will  modify 
equation  (2.40)  to  the  following: 


AA(0)U 


:(°) 


{A  Eln 


(2.66) 


Where  the  superscript  0  indicates  the  baseline  configuration.  Using  the 
process  described  for  ABC  FRF  development,  we  recalculate  the  baseline 
damaged  and  baseline  undamaged  structure  to  include  an  ABC: 


{A£/„ 


(2.67) 


Where  the  superscript  (1)  indicates  a  specific  ABC  condition  like  a  pin 
connection  on  the  first  node.  This  equation  is  then  appended  onto  the  bottom  of 
equation  (2.66) 


(2.68) 


This  process  is  then  repeated  to  the  damaged  and  undamaged  beams  as 
many  times  as  necessary  to  develop  a  full  rank,  square  matrix. 


{at(0)} 

.  K’} 

{a/L(,,)} 

Ultimately,  because  the  damage  and  change  in  frequencies  are  known, 
we  can  input  the  new  resonance  frequencies  and  sensitivity  matrices,  as 
determined  in  equation  (2.68),  into  equation  (2.46)  to  prove  that  the  new 
sensitivity  values  can  find  damage  in  all  elements. 

The  thesis  presented  by  Damanakis  [1 1 , 24]  took  this  a  step  further.  For  a 
specific  beam,  he  was  able  to  utilize  the  QR  decomposition  of  each  sensitivity 
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matrix  and  extract  specific  rows  from  a  large  variety  of  ABC  conditions  that 
provided  a  good  correlation  for  damage  detection  for  a  specific  element.  This 
method  also  develops  a  full  rank,  square  matrix  and  has  proven  reliable  in 
damage  detection. 

For  this  thesis,  we  are  exploring  the  placement  of  a  massless  spring  of 
variable  stiffness  or  a  lumped  mass  of  variable  value  on  a  single  node  of  the 
damaged  and  undamaged  beam.  Then,  using  the  process  described  in  this 
chapter,  develop  a  new  composite  sensitivity  matrix  for  damage  detection  and 
model  updating. 

By  placing  a  spring  or  mass  on  the  node  as  opposed  to  a  pin  connection, 
we  can  provide  a  near  infinite  amount  of  new  sensitivity  matrices.  However,  this 
new  method  needs  to  be  verified  that  the  change  in  resonance  frequencies  is 
sufficient  enough  to  provide  a  well  conditioned  sensitivity  matrix  that  will  detect 
damage  accurately. 

G.  CONCLUDING  REMARKS 

As  has  been  detailed  in  this  chapter,  this  thesis  will  be  using  a  2- 
dimensional  Euler-Bernoulli  beam.  In  an  effort  to  detect  damage  or  update  the 
FEM,  a  massless  spring  of  variable  stiffness  or  a  lumped  mass  of  variable  value 
has  been  introduced  into  the  system  to  provide  new  resonance  frequencies  and 
sensitivity  matrices  based  on  the  theory  of  Synthesis  and  the  Eigenvector 
Sensitivity  method. 


32 


III.  EXPERIMENTAL  SETUP  AND  TESTING 


A.  EXPERIMENTAL  SETUP 

The  structure  used  in  this  thesis  was  a  6061  Aluminum  beam  with  the 
properties  detailed  in  Table  1  and  Figure  8: 


Table  1 .  Properties  of  6061  Aluminum  Beam 


Length  cm  (in) 

182.88  (72) 

Width  cm  (in) 

4.9  (1.93) 

Height  cm  (in) 

0.947  (0.373) 

Modulus  of  elasticity  MPa 
(psi) 

68947.5729 

(10000000) 

Node/Accerometer 

Locations  mm(in) 

1.27(0.5),  14.15(5.57), 

27.03(10.64),  39.90(15.71), 

52.78(20.78),  65.66(25.85), 

78.54(30.92),  91.41(35.99), 

104.29(41.06),  117.17(46.13), 
130.05(51.20),  142.93(56.27), 
155.80(61.34),  168.68(66.41), 
181.56(71.48) 

1  2  3  4  5  6  7  8  9  10  11  12  13  14  15 


Figure  8.  Experimental  Beam  Node  Numbering 

To  reduce  complications  with  joint  modeling  and  to  be  consistent  with 
theory  developed  in  Chapter  II,  the  beam  was  set  up  in  a  free-free  condition  by 
suspending  the  beam  from  the  ceiling.  To  ensure  the  pendulum  frequency  would 
not  be  present  in  the  harmonic  response  of  the  beam,  the  beam  was  suspended 
with  an  elasto-polymer  cord  at  a  length  of  approximately  270  cm  (106.3  in).  At 
this  length,  the  idealized  pendulum  frequency  (0«O)  is  approximately  0.54Flz. 
Additionally,  the  elasto-polymer  cord  is  made  of  a  soft  core  material  that  is 
significantly  less  weight  when  compared  to  the  beam  to  accommodate  the  heave 
rigid  body  motion  of  the  cord. 
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The  beam  response  was  measured  using  Pulse  Reflex  software,  version 
21,  as  developed  by  the  Bruel  Kjaer  company  [Pulse].  Instrumentation  secured 
onto  the  beam  was  a  Piezotronics,  Inc.,  model  336C04  ICP  Accelerometer  with 
the  properties  detailed  in  Table  2  and  Table  3: 


Table  2.  Accelerometer  Properties 


Screw  Weight  g  (lb) 

0.5  (0.0011  ) 

Accelerometer  weight  g  (lb) 

6.5  (0.0143) 

Accelerometer  height  cm  (in) 

1.236  (0.487) 

Accelerometer  depth  cm  (in) 

1.574  (0.620) 

Table  3.  Accelerometer  Locations 


Node 

Serial 

Number 

Calibration 

Date 

Voltage 

sensitivity  (mV/ 
9) 

1 

10867 

9/1 0/1 993 

98.4 

2 

10875 

9/1 0/1 993 

99.3 

3 

10857 

9/1 0/1 993 

99.3 

4 

10851 

9/1 0/1 993 

98.9 

5 

10858 

9/1 0/1 993 

98.4 

6 

10859 

9/1 0/1 993 

98.7 

7 

10856 

9/1 0/1 993 

98.8 

8 

10854 

9/1 0/1 993 

95.5 

9 

11798 

2/1 7/1 993 

97.1 

10 

10877 

9/1 0/1 993 

96.4 

11 

10226 

7/27/1 993 

99.2 

12 

10847 

9/1 0/1 993 

97.9 

13 

10860 

9/1 0/1 993 

99.1 

14 

10868 

9/1 0/1 993 

96.5 

15 

10866 

9/1 0/1 993 

99.0 

Piezotronics,  Inc.  also  manufactured  the  load  cells  used  in  testing  with  the 
specifications  detailed  in  Table  4: 
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Table  4.  Load  Cell  Specifications 


MN 

SN 

Cal 

mV/g  (mV/lb) 

Hammer 

086B03 

2195 

5/22/1986 

0.02089  (9.48) 

Shaker 

208A02 

10875 

2/16/1994 

0.17809 

(80.78) 

Because  we  are  collecting  data  specifically  for  the  location  of  resonance 
peaks,  the  current  calibration  dates  of  the  accelerometers  and  load  cells  are 
acceptable  as  long  as  measurements  are  repeatable.  A  more  current  calibration 
would  be  required  if  a  study  was  being  conducted  on  waveforms  or  effects  on 
amplitudes  of  resonance  peaks. 

As  identified  in  Figure  8,  the  Accelerometers  are  offset  by  12.7mm  (0.5in) 
from  the  ends  of  the  beam.  This  was  done  to  ensure  the  accelerometer  remained 
flush  during  the  measurement  and  is  accounted  for  in  Chapter  IV 

B.  EXPERIMENTAL  TESTING 

The  beam  was  excited  using  four  different  modes:  Single  Impact  Hammer, 
Random  Impact  Hammer,  Random  Burst  Shaker,  Sine-Sweep  Shaker 

The  single  impact  hammer  measurement  utilizes  a  calibrated  impact 
hammer  to  excite  the  beam  at  a  node  and  accelerometers  measure  the  response 
over  a  given  length  of  time  with  a  single  strike.  This  process  is  repeated  several 
times  per  node  and  averaged,  and  repeated  for  each  node. 

The  Random  Impact  Hammer  test  is  conducted  similarly  to  the  single 
impact  hammer,  but  has  multiple  strikes  taken  over  time,  providing  a  larger 
amount  of  excitation.  Like  the  single  impact  hammer  test,  the  random  impact 
hammer  measurement  is  taken  for  each  node  to  formulate  a  full  FRF  matrix 

The  Burst  Random  Shaker  test  utilizes  an  electrodynamic  shaker  to  excite 
the  beam  with  a  randomly  generated  burst  of  frequencies  and  measure  the 
response.  This  test  requires  a  stand  to  isolate  the  shaker  from  outside 

frequencies  while  it  conducting  the  measurement  and  to  be  attached  to  the  beam 
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with  a  sufficiently  stiff  “stinger”  to  reduce  the  loss  of  energy  generated  by  the 
shaker.  Similar  to  the  hammer  tests,  this  process  is  repeated  for  each  node  to 
formulate  a  full  FRF  matrix. 

The  Sine-Sweep  shaker  test  also  utilizes  an  electrodynamic  shaker  to 
excite  the  beam.  However,  instead  of  exciting  the  beam  with  a  very  short  burst  of 
random  frequencies,  the  sine-sweep  method  excited  the  beam  at  every 
frequency  given  a  specified  bandwidth  and  step.  This  is  done  slowly  as  to  ensure 
the  beam  is  responding  only  to  the  given  input  forcing  function  and  the 
accelerometers  can  accurately  measure  the  response.  Again,  this  process  is 
repeated  for  each  node  to  formulate  a  full  FRF  matrix 

Ultimately,  the  single  impact  hammer  test  and  random  impact  hammer  test 
proved  far  more  accurate  and  repeatable  than  the  burst  random  shaker  test  and 
the  sine-sweep  shaker  test.  The  shaker  tests  where  not  useable  due  to  the  mass 
ratio  of  the  shaker  to  the  beam  significantly  altering  the  location  of  the  resonance 
frequencies  as  displayed  in  Figure  9. 
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Random  Impact 
Random  Burst 
Single  Impact 
Step  Sine 


Graph  developed  using  [7] 

Figure  9.  Testing  Signal  Comparison 


Using  the  measurements  from  the  single  impact  hammer,  the 
measurement  software  [7]  was  able  to  determine  the  resonance  frequencies  of 
the  beam  using  a  rational  fraction  polynomial,  parameter  estimation  method.  The 
results  of  which  are  displayed  in  Figure  10. 


37 


Mode  Selection 


Graph  created  using  [7], 

Figure  10.  Single  Impact  Hammer  Testing  Mode  Selection 

C.  CONCLUDING  REMARKS 

Experimental  testing  using  the  impact  hammer  provided  much  clearer 
results  in  determining  the  modal  response  of  the  beam.  Tests  using  the  shaker 
setup  did  not  provide  a  clear  picture  of  the  lower  frequencies  possibly  due  to  the 
mass  ratio.  Investigations  and  research  should  be  conducted  to  determine  the 
limits  of  using  the  electro-dynamic  shaker  because  it  is  capable  of  producing  high 
quality  results  with  a  dramatic  decrease  in  time  given  the  correct  setup. 
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IV.  FINITE  ELEMENT  MODELING  AND  SIMULATION 


A.  FINITE  ELEMENT  MODELING 

The  Finite  Element  Model  was  constructed  using  16,  2-D  Euler-Bernoulli 
elements.  The  FEM  had  17  nodes  placed  on  either  end  of  the  bar  and  correlated 
to  the  placement  of  the  accelerometers  on  the  prototype  beam.  As  such,  the 
FEM  had  the  properties  displayed  in  Table  5  and  Figure  1 1 : 


Table  5.  Properties  of  Beam  FEM 


Length  cm  (in) 

182.88  (72) 

Width  cm  (in) 

4.9  (1.93) 

Height  cm  (in) 

0.947  (0.373) 

Modulus  of  elasticity  MPa 

68947.5729 

(psi) 

(10000000) 

Density  Kg/m'3  (Ibm/in3) 

2698.79  (0.0975) 

Node  Locations  cm(in) 

0(0),  1.27(0.5),  14.15(5.57), 
27.03(10.64),  39.90(15.71), 
52.78(20.78),  65.66(25.85), 
78.54(30.92),  91.41(35.99), 
104.29(41.06),  117.17(46.13), 
130.05(51.20),  142.93(56.27), 
155.80(61.34),  168.68(66.41), 
181.56(71.48),  182.83(71.98) 

Figure  1 1 .  FEM  Node  Numbering 


In  total,  the  FEM  has  34  DOF,  where  odd  numbered  DOF  are  related  to 
translation  and  even  numbered  DOF  are  related  to  rotation  as  identified  in  Figure 
12. 
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Figure  12.  2-D  Euler-Bernouili  Element  DOF  Setup.  Adapted  from  [25]. 


In  addition  to  formulating  the  Mass  and  Stiffness  properties  for  the  beam 
with  the  data  from  Table  5,  the  instrumentation  for  the  accelerometers  where 
added  to  the  mass  matrix  as  a  point  load  with  inertial  mass.  As  such,  the  nodes 
that  are  not  on  the  ends  have  the  instrumentation  loads  added  to  the  diagonal. 

In  addition  to  the  FEM,  a  second  model  was  constructed  using  the  same 
specifications  in  Table  2  and  Table  5,  but  damaged.  The  damage  was  simulated 
by  reducing  the  flexural  rigidity  (El)  in  a  single  element. 

Additional  aspects  for  setting  up  the  FEM  are  detailed  in  Appendix  A. 

B.  FINITE  ELEMENT  MODEL  TESTING 

Following  the  theory  in  Chapter  II,  the  FEM  mass  and  stiffness  matrices 
are  constructed  and  used  to  solve  the  eigenvalue  problem  as  described  for 
equation  (2.18).  Using  the  “eig”  function  in  MATLAB  [17],  we  calculate  a  series  of 
resonance  frequencies  and  mass-normalized  mode  shapes.  To  confirm  their 
accuracy,  the  frequencies  are  compared  to  an  analytic  solution  using  the  roots  to 
the  characteristic  equation  for  a  uniform,  cantilevered  beam  utilizing  equation 
(4.1)  from  [4]  : 


(III  _ 

L  JpA  ’ 


(4.1) 
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where  (XL,)  is  a  tabulated  value  found  in  [26]  for  specific  modes  given  specific 

boundary  conditions.  This  thesis  used  the  following  values  to  predict  the  first  4 
modes: 

TLj  =1.8751 
TL,  =  4.6941 

(4.2) 

/IL3  =  7.8548 
AL4=  10.996 

Table  6  displays  the  accuracy  of  the  FEM  to  an  analytical  assessment: 


Table  6.  Analytical  and  FEM  Frequency  Comparison 


Analytical 

Frequency 

FEM 

Frequency 

%  Rel 

Error 

Mode  1 

2.22641373 

2.2276613 

0.00163 

Mode  2 

13.9528099 

13.960624 

0.01021 

Mode  3 

39.0685476 

39.092269 

0.03098 

Mode  4 

76.5643263 

76.618568 

0.07084 

Because  the  percent  relative  error  is  less  than  0.1%,  we  can  say  the  FEM 
has  good  correlation  to  analytical  values.  The  FEM  is  capable  of  calculating  34 
total  natural  frequencies,  however,  the  equipment  measuring  the  response  of  the 
beam  is  limited  to  1000  Hz,  which  captures  the  first  12  modes  of  the  beam. 
Additionally,  because  we  are  using  a  beam  in  a  free-free  boundary  condition,  the 
first  two  calculated  modes  are  considered  rigid  body  modes  due  to  their  proximity 
to  zero  Hz  resonance  and  their  mode  shape.  For  the  purposes  of  calculations, 
rigid  body  modes  are  removed.  Table  7  compares  the  FEM  resonance  values  to 
the  measured  values: 
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Table  7.  FEM  Resonance  Frequencies 


FEM 

Frequency  (Hz) 

Measured 

Frequency 

(Hz) 

Absolute 
Error  (Hz) 

%  Rel 

Error 

Mode  1 

14.35143 

13.81281 

0.53862 

3.8994 

Mode  2 

39.55544 

38.27358 

1.28186 

3.3492 

Mode  3 

77.54409 

75.02992 

2.51417 

3.3509 

Mode  4 

128.20934 

124.01327 

4.19607 

3.3836 

Mode  5 

191.61915 

185.21543 

6.40372 

3.4574 

Mode  6 

267.88319 

258.58254 

9.30065 

3.5968 

Mode  7 

357.17924 

344.25651 

12.92273 

3.7538 

Mode  8 

459.76188 

442.25887 

17.50301 

3.9576 

Mode  9 

575.95281 

552.67146 

23.28135 

4.2125 

Mode  10 

706.07569 

674.52174 

31.55395 

4.678 

Mode  11 

850.20479 

808.39895 

41.80584 

5.1714 

Mode  12 

1007.10225 

953.93877 

53.16348 

5.573 

To  adjust  for  the  discrepancies,  a  3%  correction  factor  used  on  the 
modulus  of  elasticity.  The  results  of  which  are  displayed  in  Table  8. 


Table  8.  FEM  Resonance  Frequencies  with  Correction  Factor 


Corrected  FEM 
Frequency(Hz) 

Measured 
Frequency  (Hz) 

Absolute 
Error  (Hz) 

%  Rel 

Error 

Mode  1 

14.13452 

13.81281 

0.32171 

2.329 

Mode  2 

38.95759 

38.27358 

0.68401 

1.7872 

Mode  3 

76.37207 

75.02992 

1.34215 

1.7888 

Mode  4 

126.27156 

124.01327 

2.25829 

1.821 

Mode  5 

188.72298 

185.21543 

3.50755 

1.8938 

Mode  6 

263.83435 

258.58254 

5.25181 

2.031 

Mode  7 

351.78075 

344.25651 

7.52424 

2.1856 

Mode  8 

452.81294 

442.25887 

10.55407 

2.3864 

Mode  9 

567.24773 

552.67146 

14.57627 

2.6374 

Mode  10 

695.40391 

674.52174 

20.88217 

3.0958 

Mode  11 

837.35461 

808.39895 

28.95566 

3.5819 

Mode  12 

991.88069 

953.93877 

37.94192 

3.9774 
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The  new  FEM  values  correlate  to  previous  work  [11]  and  prove 
consistency  among  the  models. 

C.  FEM  INITIAL  SENSITIVITY  MATRIX  AND  FRF 
1.  Initial  Sensitivity  Matrix 

The  initial  sensitivity  matrix  is  influence  by  the  number  of  modes  that  are 
available  for  calculation.  For  this  FEM,  if  we  use  all  modes  that  can  be 
calculated,  then  we  will  get  an  over  determined  sensitivity  matrix  of  34  by  14. 
However,  because  of  the  equipment  bandwidth  limitation  of  1,000  hertz,  we  only 
use  the  first  12  flexural  modes  of  the  initial  matrix  for  calculation.  As  such,  the 
initial  sensitivity  matrix  is: 


2.33S-05 

0.0004 

0.0021 

0.0054 

0.0098 

0.0141 

0.0168 

0.0168 

0.0141 

0.0098 

0.0054 

0.0021 

0.0004 

0.0000 

0.0011 

0.0178 

0.0631 

0.1083 

0.1088 

0.0604 

0.0105 

0.0105 

0.0604 

0.1088 

0.1083 

0.0631 

0.0178 

0.0011 

0.0144 

0.1751 

0.4148 

0.3483 

0.0709 

0.0743 

0.3244 

0.3244 

0.0743 

0.0709 

0.3483 

0.4148 

0.1751 

0.0144 

0.0900 

0.8146 

1.0697 

0.2027 

0.4223 

0.9942 

0.2939 

0.2939 

0.9942 

0.4223 

0.2027 

1.0697 

0.8146 

0.0900 

0.3710 

2.3396 

1.2645 

0.6437 

2.1226 

0.3942 

1.5482 

1.5482 

0.3942 

2.1226 

0.6437 

1.2645 

2.3396 

0.3710 

1.1525 

4.7002 

0.8986 

3.6816 

1.3456 

3.0337 

2.1594 

2.1594 

3.0337 

1.3456 

3.6816 

0.8986 

4.7002 

1.1525 

2.9140 

7.0920 

2.4919 

5.5831 

3.2166 

4.8527 

4.0214 

4.0214 

4.8527 

3.2166 

5.5831 

2.4919 

7.0920 

2.9140 

6.2837 

8.6226 

8.2602 

4.0817 

10.4211 

4.2973 

8.0243 

8.0243 

4.2973 

10.4211 

4.0817 

8.2602 

8.6226 

6.2837 

11.9021 

9.7756 

14.9724 

8.9028 

8.8994 

15.1000 

8.8989 

8.8989 

15.1000 

8.8994 

8.9028 

14.9724 

9.7756 

11.9021 

20.1484 

13.1424 

17.9452 

20.1283 

14.1650 

13.1212 

19.2531 

19.2531 

13.1212 

14.1650 

20.1283 

17.9452 

13.1424 

20.1484 

30.6023 

21.6249 

21.7209 

24.8308 

26.3369 

24.3565 

21.4784 

21.4784 

24.3565 

26.3369 

24.8308 

21.7209 

21.6249 

30.6023 

40.5225 

33.6458 

34.8236 

34.8729 

33.6127 

31.8271 

30.5621 

30.5621 

31.8271 

33.6127 

34.8729 

34.8236 

33.6458 

40.5225 

As  detailed  for  equation  (2.31),  the  sensitivity  matrix  is  sized  for  the 
number  of  modes  by  the  number  of  interested  elements.  Because  we  are  limited 
to  12  modes,  the  initial  sensitivity  matrix  is  sized  14  elements  by  12  modes. 
Additionally,  because  this  sensitivity  matrix  is  developed  for  a  free-free  boundary 
condition,  the  matrix  is  symmetric  and  thus  only  rank  6.  The  under-determined, 
rank  deficient,  initial  matrix  is  not  suitable  for  damage  detection  or  model 
updating  and  will  need  additional  modes  and  information  before  it  can  be  used  as 
such. 


As  a  proof,  equation  (2.32)  is  used  to  check  the  sensitivity  matrix  for 
accuracy  by  using  the  sensitivity  value  to  predict  a  new  resonance  frequency  for 
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a  given  change  in  El  at  a  specified  element,  and  compared  the  resonance 
frequency  calculated  from  an  FEM  with  the  change  installed. 


‘0.0002 

0.0043 

0.0169 

0.0341 

0.0463 

0.0450 

0.0283 

0.0014 

0.0256 

0.0428 

0.0448 

0.0334 

0.0166 

0.0042 

0.0016 

0.0225 

0.0614 

0.0618 

0.0019 

0.0643 

0.0666 

0.0001 

0.0674 

0.0661 

0.0004 

0.0598 

0.0603 

0.0223 

0.0051 

0.0566 

0.0848 

0.0222 

0.0966 

0.0014 

0.0881 

0.0008 

0.0870 

0.0010 

0.0977 

0.0242 

0.0830 

0.0559 

0.0116 

0.0933 

0.0338 

0.1105 

0.0284 

0.0740 

0.0891 

0.0003 

0.0903 

0.0726 

0.0278 

0.1117 

0.0320 

0.0921 

0.0214 

0.1136 

0.0608 

0.0353 

0.0854 

0.0987 

0.0666 

0.0006 

0.0658 

0.0997 

0.0843 

0.0360 

0.0622 

0.1122 

0.0340 

0.1048 

0.1111 

0.0822 

0.0681 

0.0500 

0.0252 

0.0004 

0.0261 

0.0491 

0.0691 

0.0812 

0.1122 

0.1036 

0.0483 

0.0697 

0.0757 

0.0516 

0.0388 

0.0275 

0.0134 

0.0004 

0.0142 

0.0267 

0.0396 

0.0508 

0.0766 

0.0687 

0.0628 

0.0240 

0.0033 

0.0412 

0.0637 

0.0607 

0.0376 

0.0004 

0.0368 

0.0616 

0.0629 

0.0421 

0.0041 

0.0231 

0.0759 

0.0128 

0.0335 

0.0383 

0.0004 

0.0399 

0.0391 

0.0004 

0.0399 

0.0391 

0.0004 

0.0391 

0.0327 

0.0139 

0.0854 

0.0290 

0.0208 

0.0095 

0.0249 

0.0039 

0.0265 

0.0003 

0.0257 

0.0049 

0.0257 

0.0089 

0.0201 

0.0302 

0.0895 

0.0257 

0.0009 

0.0096 

0.0047 

0.0056 

0.0081 

0.0006 

0.0090 

0.0061 

0.0042 

0.0088 

0.0002 

0.0268 

0.0845 

0.0139 

0.0029 

0.0007 

0.0020 

0.0032 

0.0023 

0.0002 

0.0028 

0.0040 

0.0031 

0.0006 

0.0018 

0.0148 

The  largest  difference  is  0.1136%  relative  error  occurring  at  element  2, 
mode  5.  As  this  error  is  less  than  1%,  we  can  infer  the  sensitivity  matrix  has  good 
prediction  qualities. 

2.  Initial  Calculated  FRF 

Using  the  process  detailed  in  Chapter  II,  the  initial  FEM  FRF  is  calculated 
and  plotted  in  Figure  13.  The  initial  FEM  FRF  is  then  compared  to  the  measured 
FRF  in  Figure  14. 


Frequency (Hz) 

Figure  13.  FEM  Initial  FRF  Response,  Hi,i 
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Figure  14.  Initial  FEM  vs  Measured  FEM  Response,  Flu 


In  comparison  with  the  measured  response  as  detailed  in  Chapter  III,  the 
FEM  and  the  model  do  not  match  up  as  displayed  in  Figure  14.  This  indicates 
that  the  model  will  need  to  undergo  model  updating  to  correctly  reflect  the 
prototype  structure. 

The  proposed  method  for  model  updating  and  damage  detection  is  to 
synthesize  a  massless  spring  or  a  lumped  mass  into  the  FEM  to  develop  new 
sensitivity  matrices  that  can  be  used  to  develop  a  composite  sensitivity  matrix  for 
damage  detection  and  model  updating.  The  first  step  of  this  is  to  see  if  these 
variables  can  be  synthesized  into  the  new  system 
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D.  SINGLE  SPRING  DYNAMIC  RESPONSE 

As  stated  earlier,  one  of  the  benefits  of  synthesis  is  the  capability  to 
manipulate  data  without  the  tedious  time  and  expense  of  actually  installing  the 
boundary  conditions  on  the  structure  and  performing  another  vibration  test. 

1.  FRF  Synthesis 

Previous  work  [5,  11,  24]  utilized  a  simulated  pin  to  modify  the  structure 
response.  For  example,  following  the  theory  and  procedure  in  Chapter  II,  if  a  pin 
was  placed  on  node  13  of  the  beam,  we  would  expect  the  response  plotted  in 
Figure  15. 


Figure  15.  Initial  FRF  Response  vs  Synthesized  Pin  Located  at  Node  13, 

H33,33 

To  ensure  that  the  synthesis  is  accomplished  correctly,  a  second  FRF  is 
calculated  with  the  pin  installed.  This  removes  the  associated  row  and  column  for 
the  DOF  the  pin  is  attached  to.  The  synthesis  check  is  compared  in  Figure  16. 
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Figure  16.  Pin  FRF  Response  Check,  Pin  Located  at  Node  13,  H33;33 

However,  if  this  approach  was  unable  to  produce  enough  rows  for  a  much 
larger  structure,  damage  detection  and  model  updating  could  not  be  reliably 
accomplished.  As  such,  the  synthesis  of  a  variable  quantity  should  provide  more 
information 

2.  Single  Spring  FRF  Synthesis 

The  introduction  of  a  massless  spring  introduces  the  possibility  of  an 
infinite  amount  of  new  data  as  the  spring  can  be  any  variable,  and  can  be  placed 
at  any  node. 

To  synthesize  the  spring,  we  follow  the  procedure  detailed  in  Chapter  II 
and  identify  the  CSET  as  the  translational  DOF  for  the  node  the  spring  is  placed 
on.  Similar  to  the  pin  synthesis,  we  also  perform  a  check  of  this  synthesis  to 
ensure  an  exact  prediction  as  detailed  in  Figure  17. 
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Figure  1 7.  1 78.579  kg/mm  (1 0,000lb/in)Spring  Synthesis  at  Node  2 

Check,  H3i3 

Using  equation  (2.64)  and  (2.65),  we  are  able  to  predict  the  modal 
response  of  the  beam  with  the  spring  installed  as  detailed  in  Figure  18  and 
Figure  19: 


Spring  1  Lb 
Spring  lOLb 
Spring  lOOLb 
Spring  500Lb 
Spring  lOOOLb 
Spring  lOOOOLb 

-  Pin 

-  Inital 


Figure  18.  FRF  Progression  for  Spring  Synthesis  at  Node  13,  H25,25 
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Frequency  (Hz) 


Spring  ILb 
Spring  lOLb 
Spring  lOOLb 
Spring  500Lb 
Spring  lOOOLb 
Spring  lOOOOLb 
Pin 
Inital 


Figure  19.  Magnification  of  FRF  Progression  for  Spring  Synthesis  at  Node 

1  3,  H 25,25 


As  can  be  seen  in  Figure  19,  by  increasing  the  strength  of  the  spring,  the 
dynamic  beam  response  will  gradually  change  to  a  response  similar  to  that  of  pin 
connection.  As  such,  we  can  see  the  progression  of  the  modal  resonance 
frequencies  with  the  introduction  of  the  spring  in  Figure  20  and  Figure  21. 
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Frequency  Progression  for  Spring  Strength  0-500  Ib/in 

(Nodel3) 
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Figure  20.  Modal  Resonance  Frequency  Progression  for  Spring  Strength 
0-8.928  kg/mm  (0— 500lb/in)  Located  at  Node  13 
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Figure  21 .  Modal  Resonance  Frequency  Progression  for  Spring  Strength 
0-178.579  kg/mm  (0—1 0,0001b/in)  Located  at  Node  13 


50 


Even  at  low  spring  strengths,  we  see  a  loss  of  one  of  the  two  rigid  body 
modes.  To  understand  what  the  beam  response  is  at  these  points,  the  mode 
shapes  where  mapped  out  in  Figure  22  through  Figure  31 : 


0.02 


— Model(O)  — Mode2(0)  —•—Mode  3(13.86)  —•—Mode  4  (38.21)  A  Sprir^  Location 


Figure  22.  Mode  Shapes  with  0  kg/mm(0  Ib/in)  Spring  Located  at  Node  1 3 

(Resonance  Frequency  in  Parentheses) 


0.02 


■Mode  1(0) 


Figure  23. 


— ♦-  Mode  2  (1.92)  Mode  3  (13.86)  Mode  4  (38.23)  A  Spring  Location 

Mode  Shapes  with  0.017  kg/mm  (1  Ib/in)  Spring  Located  at 
Node  13  (Resonance  Frequency  in  Parentheses) 
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Node  Vertical  Position  -n  Node  Vertical  Position 


0.02 


•Mode  1(0)' 


•Mode  2  (8.46) 


■Mode  3  (13.87)' 


•  Mode  4  (38.68)  A  Spring  Location 


:igure  24.  Mode  Shapes  with  0.357  kg/mm(20  Ib/in)  Spring  Located  at 
Node  13  (Resonance  Frequency  in  Parentheses) 


0.02 


Mode  1  (0)  — Mode  2  (12.92)  — Mode  3  (13.91)  — Mode  4  (39.42)  A  Spring  Location 

Figure  25.  Mode  Shapes  with  0.892  kg/mm(50  ibf/in)Spring  Located  at 
Node  13  (Resonance  Frequency  in  Parentheses) 
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Node  Vertical  Position  -n  Node  Vertical  Position 


0.02 

0.015 

0.01 

0.005 

0 

-0.005 

-0.01 

-0.015 


-0.02 


-0.025 


-0.03 
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■  Mode  1  (0)  ■ 


•Mode  2  (13.42)  • 


■  Mode  3  (13.97)  ■ 


■  Mode  4  (39.55)  A  Spring  Location 


:igure  26.  Mode  Shapes  with  0.982  kg/mm(55  Ib/in)  Spring  Located  at 
Node  13  (Resonance  Frequency  in  Parentheses) 


0.02 


— •—  Mode  1  (0)  — Mode  2  (13.50)  — •—  Mode  3  (13.99)  Mode  4  (39.57)  A  Spring  Location 


Figure  27.  Mode  Shapes  with  1 .00  kg/mm  (56  Ib/in)  Spring  Located  at 
Node  13  (Resonance  Frequency  in  Parentheses) 
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Node  Vertical  Position  -n  Node  Vertical  Position 
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•  Mode  1  (0)  ■ 


•Mode  2  (13.84) 


•  Mode  3  (17.40)  • 


-Mode 4  (40.73)  A  Spring  Location 


igure  28.  Mode  Shapes  with  1 .79  kg/mm  (100  ib/in)  Spring  Located  at 
Node  13  (Resonance  Frequency  in  Parentheses) 


0.02 


— •—  Mode  1  (0)  — Mode  2  (13.85)  — •—  Mode  3  (26.69)  Mode  4  (49.51)  A  Spring  Location 


Figure  29.  Mode  Shapes  with  8.93  kg/mm(500  lbf/in)Spring  Located  at 
Node  13  (Resonance  Frequency  in  Parentheses) 
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0.02 


-♦-Mode  1  (0)  — ♦—  Mode  2  (13.85)  — ♦—  Mode  3  (26.69)  — ♦—  Mode  4  (49.51)  A  Spring  Location 

Figure  30.  Mode  Shapes  with  178.579  kg/mm  (10,000  Ib/in)  Spring 
Located  at  Node  13  (Resonance  Frequency  in  Parentheses) 


0.02 


—♦—Mode  1(0)  —♦—Mode  2  (13.85)  —♦—Mode  3  (30.15) —♦—Mode  4  (59.17)  A  Pin  Location 


Figure  31 .  Mode  Shapes  with  a  Pin  Connection  Located  at  Node  1 3 
(Resonance  Frequency  in  Parentheses) 

As  Figure  22  through  Figure  31  display,  by  adding  a  spring  to  the  beam, 
even  at  small  values,  we  transfer  a  rigid  body  mode  into  a  flexural  mode.  As  the 
spring  increases  in  strength,  the  mode  shapes  gradually  transform  the  next  lower 
mode  (mode  4  transforms  to  mode  3).  As  the  spring  approaches  a  near  infinite 
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strength,  the  dynamic  beam  response  resembles  that  of  an  artificial  pin 
connection  placed  at  the  selected  node.  To  see  if  this  transformation  ultimately 
provides  new  and  relevant  information  for  damage  detection  and  model  updating 
we  looked  at  the  values  of  the  sensitivity  matrix. 

3.  Single  Spring  Sensitivity  Analysis 

As  discussed  in  Chapter  II,  the  sensitivity  matrix  is  what  is  ultimately  used 
in  damage  detection  and  model  updating.  Previous  thesis  work  [11]  indicated 
higher  sensitivity  modes  provide  greater  fidelity  in  damage  detection.  As  such, 
we  tracked  the  beam  sensitivity  values  of  the  highest  mode  in  the  sensitivity 
matrix  as  the  spring  increased  from  0-10,000  Ibf/in. 


[Mode(l)] 
| ~  Mode  (2)] 


[Modern-  l)j 
| 'Mode  (ft)] 


<-  HighestMode 


(4.5) 


It  is  important  to  note  that,  for  a  majority  of  these  simulations,  the 
resonance  frequency  for  the  highest  mode  goes  above  the  1,000  hertz  limit  at 
some  point.  When  this  happens,  the  sensitivity  matrix  no  longer  calculates  the 
sensitivity  to  that  mode,  so  the  highest  mode  shifts  to  the  next  lower  frequency. 
Ultimately,  this  will  assist  in  damage  detection  and  model  updating.  The  results  of 
the  simulations  are  displayed  in  Figures  32  though  Figure  46. 
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igure  32. 


Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  2 
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Figure  33.  Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  3 
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igure  34. 


Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  4 
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Figure  35.  Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  5 
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Figure  36.  Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  6 
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Figure  37.  Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  7 
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:igure  38. 


Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  8 
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Figure  39.  Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  9 
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:igure  40. 


Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  10 
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Figure  41 .  Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  1 1 
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Figure  42.  Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  12 
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Figure  43.  Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  13 


62 


50 


5 


0 

0.0  20.0  40.0  60.0  80.0  100.0  120.0  140.0  160.0  180.0 

Spring  Strength  (kg/mm) 

Figure  44.  Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  14 
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Figure  45.  Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  15 
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Figure  46.  Effects  of  Spring  Strength  Variance  on  Sensitivity  Value  for 

Node  16 

As  was  alluded  to  in  Chapter  II,  one  of  the  concerns  of  this  method  is 
matrix  conditioning,  otherwise  known  as  the  matrix  condition  number.  For  the 
purposes  of  this  thesis,  the  condition  number  relates  to  the  linear  dependence  of 
the  rows  and  columns  in  the  matrix.  An  example  of  how  this  effects  our  capability 
to  locate  damage  is  as  follows: 

If  we  use  the  all  of  the  baseline  sensitivity  matrix  and  supplement  two 
additional  modes  from  a  sensitivity  matrix  that  has  an  8.928  kg/mm  (500lb/in) 
spring  installed  at  node  9,  we  get  a  square,  full  rank  matrix  that  has  a  matrix 
condition  number  of  6.7441  e+11.  This  high  of  a  condition  number  indicates  a 
poorly  conditioned  matrix  that  most  likely  will  not  produce  accurate  results.  For 
example,  Figure  47  is  the  result  of  using  this  sensitivity  matrix  to  find  5%  damage 
at  node  6 
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Figure  47.  Model  Error  vs  True  Error  Using  Poorly  Conditioned  Sensitivity 

Matrix 


However,  sensitivity  matrix  that  has  two  baseline  modes  and  twelve 
modes  from  a  sensitivity  matrix  that  has  a  178.579  kg/mm  (1  OOOOIb/in)  spring 
installed  at  node  2,  we  also  have  a  square,  full  rank  matrix  with  a  matrix  condition 
number  of  1.7671e+04.  This  condition  number  indicates  a  better  conditioned 
matrix  that  most  likely  will  produce  accurate  results.  For  example,  Figure  48  is 
the  result  of  using  this  matrix  to  find  the  same  5%  damage  at  node  6. 
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Figure  48.  Model  Error  vs  True  Error  Using  a  Well  Conditioned  Sensitivity 

Matrix 


Because  of  this,  we  identified  spring  stiffness  values  from  Figures  32 
through  Figure  46  that  varied  greatly  from  baseline,  and  from  each  other.  These 
spring  strengths  are  identified  in  Table  9. 
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Table  9.  Useable  Single  Spring  Strengths 


Node 

Spring  Strength  kg/mm  (Ib/mm) 

2 

71.43(4000) 

178.58(10000) 

3 

107.15(6000) 

178.58(10000) 

4 

64.29(3600) 

178.58(10000) 

5 

71.43(4000) 

178.58(10000) 

6 

67.86(3800) 

178.58(10000) 

7 

116.08(6500) 

178.58(10000) 

8 

71.43(4000) 

178.58(10000) 

9 

Not  useable 

10 

71.43(4000) 

178.58(10000) 

11 

116.08(6500) 

178.58(10000) 

12 

67.86(3800) 

178.58(10000) 

13 

71.43(4000) 

178.58(10000) 

14 

64.29(3600) 

178.58(10000) 

15 

107.15(6000) 

178.58(10000) 

16 

71.43(4000) 

178.58(10000) 

As  noted  in  Table  9,  Node  9  was  determined  to  be  unusable  due  to  no 
variance  in  the  sensitivity.  This  lack  of  variance  may  be  due  to  its  location  in  the 
middle  of  the  beam  and  the  symmetric  characteristics  a  free-free  beam  exhibits. 

Because  the  baseline  sensitivity  matrix  is  both  under-determined  and  not 
full  rank,  we  used  the  sensitivity  matrices  from  the  baseline  and  those  developed 
from  the  spring  strengths  detailed  in  Table  9  to  build  an  overall  composite 
sensitivity  matrix. 


comp 


£#13-4,000 

<-,#13-10,000 


] 

] 


(4.6) 


Equation  (4.6)  is  an  example  of  how  a  composite  sensitivity  matrix  could 
be  developed  where  [s°]is  the  baseline  sensitivity  matrix  and  [5'7V13_10’000^]  is  the 
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sensitivity  matrix  as  developed  from  using  a  10,000  Ibf/in  spring  on  node  13.  To 
make  a  square  composite  matrix,  we  take  a  selected  number  of  rows  separate 
from  each  configurations.  An  easily  repeatable  method  is  to  extract  a  fixed  set 
series  of  rows  from  either  the  top  or  bottom  of  a  sensitivity  matrix  and 
determining  the  overall  matrix  condition  number.  A  high  condition  number 
indicates  poor  matrix  conditioning  and  likely  will  have  large  errors  when  used  in 
calculations,  whereas  a  low  condition  number  indicates  linear  independence 
between  the  rows  and  columns  of  the  matrix.  A  matrix  with  a  condition  number  of 
1  will  be  perfectly  accurate  in  solving  a  system  of  equations. 

As  there  are  many  different  possible  combinations  of  the  sensitivity 
matrices  that  would  still  get  a  square,  full-rank  matrix  for  damage  detection,  we 
calculated  the  condition  number  for  several  different  scenarios.  The  first  series  of 
calculations  determine  how  many  baseline  modes  should  be  used  in  the 
development  of  the  composite  sensitivity  matrix,  and  whether  to  use  higher 
frequency  modes,  or  to  use  lower  frequency  modes.  The  composite  matrix  was 
built  under  the  following  system  where  the  variable  “r”  indicates  the  number  of 
rows  from  either  the  top  (for  low  modes)  or  bottom  (for  high  modes)  of  that 
matrix.  Using  equation  (4.7),  Figure  49  and  Figure  50  show  the  effects  the 
number  of  baseline  modes  has  on  the  sensitivity  matrix  condition  number. 


L  comPJuxU 


S° 

Jrxl4 

iAf7-6,500 
» 

J(14-r)jcl4 


(4.7) 
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Figure  49.  Effect  of  Number  of  Baseline  Modes  Used  in  Sensitivity  Matrix 
on  Condition  Number,  Single  Node,  178.579  kg/mm  (10,000  lb/ 

in)  Artificial  Spring 


Figure  50.  Effect  of  Number  of  Baseline  Modes  Used  in  Sensitivity  Matrix 
on  Condition  Number,  71. 431  kg/mm  (4000lb/in)  (Node  13) 

1 16.076  kg/mm  (6500lb/in)  (Node  7) 
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For  the  development  of  Figure  51,  the  sensitivity  matrix  comprised  of 
values  from  three  matrices:  Baseline,  Single-Spring  Mid-Strength,  and  Single- 
Spring  High  Strength.  Using  the  same  rules  developing  equation  (4.7),  the 
sensitivity  matrix  was  developed  using  equation  (4.8) 


re  n 

I  comp  I 


14x14 


(4.8) 


Figure  51 .  Effect  of  Number  of  Baseline  Modes  Used  in  Sensitivity  Matrix 
on  Condition  Number,  71.431  &  178.579kg/mm 
(4,000&10,000lb/in)  on  Node  13,  116.076  &  178.579kg/mm 
(6,500&10,000lb/in)  on  Node  7 
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Figure  49  through  Figure  51  indicate  to  us  that  the  best-conditioned 
composite  sensitivity  matrices  (lowest  matrix  condition  number)  are  found  when 
we  use  two  or  less  modes  from  the  baseline  sensitivity  matrix  and  when  using  at 
least  two  modified  matrices.  Additionally,  these  figures  confirm  previous  thesis 
work  [11]  in  that  the  use  of  higher  order  modes  will  provide  a  more  accurate 
solution  in  damage  detection  as  opposed  to  using  low  order  modes. 

Using  the  identified  useable  conditions  noted  in  Table  9,  we  can  combine 
modes  from  different  spring  strengths  and  different  nodes.  We  can  combine 
these  different  conditions  into  a  composite  sensitivity  matrix,  and  then  evaluate 
the  matrix  condition  number.  From  this,  we  developed  a  map  of  possible 
combinations  for  spring  strengths  and  nodes  based  on  an  index  system,  as 
displayed  in  Table  10. 


Table  1 0.  Single  Spring  Index  for  Evaluating  Condition  Number 


Index 

Number 

Node 

Spring 

kg/mm  (Ib/in) 

1 

2 

71.43(4000) 

2 

2 

178.58(10000) 

3 

3 

107.15(6000) 

4 

3 

178.58(10000) 

5 

4 

64.29(3600) 

6 

4 

178.58(10000) 

7 

5 

71.43(4000) 

8 

5 

178.58(10000) 

9 

6 

67.86(3800) 

10 

6 

178.58(10000) 

11 

7 

116.08(6500) 

12 

7 

178.58(10000) 

13 

8 

71.43(4000) 

14 

8 

178.58(10000) 

Index 

Number 

Node 

Spring 

kg/mm  (Ib/in) 

15 

10 

71.43(4000) 

16 

10 

178.58(10000) 

17 

11 

116.08(6500) 

18 

11 

178.58(10000) 

19 

12 

67.86(3800) 

20 

12 

178.58(10000) 

21 

13 

71.43(4000) 

22 

13 

178.58(10000) 

23 

14 

64.29(3600) 

24 

14 

178.58(10000) 

25 

15 

107.15(6000) 

26 

15 

178.58(10000) 

27 

16 

71.43(4000) 

28 

16 

178.58(10000) 

71 


Using  equation  (4.8),  two  modes  where  used  from  baseline  and  six  modes 
where  used  from  two  different  indexes  identified  in  Table  10  to  form  a  square,  full 
rank  sensitivity  matrix. 
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This  composite  sensitivity  matrix  was  then  evaluated  for  its  condition 
number  and  compared  to  the  remaining  possible  combinations.  In  total,  there  are 
784  different  combinations,  all  with  different  condition  numbers  indicating  the 
accuracy  of  that  matrix’s  capability  to  evaluate  damage.  To  assist  in  determining 
what  combination  worked  best,  the  index  of  Table  10  was  used  to  develop  Figure 
52  and  Figure  53. 
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Figure  52.  Sensitivity  Matrix  Condition  Number  Index  Example 


Figure  52  is  a  magnified  view  of  all  possible  combinations  using  the  index 
provided  in  Table  10.  The  color  of  the  identified  box  indicates  the  condition 
number  of  a  composite  sensitivity  matrix  consisting  of  the  two  highest  baseline 
modes  and  of  six  highest  modes  from  the  sensitivity  matrix  from  index  1  and 
index  3  as  explained  in  equation  (4.10).  In  this  instance,  the  corresponding 
condition  number  is  approximately  2600 
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(4.10) 
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Figure  53.  Sensitivity  Matrix  Combination  Condition  Number  Indicator, 

Single  Spring 


We  want  to  minimize  the  matrix  condition  number  to  get  the  most  accurate 
damage  detection  solution.  By  analyzing  Figure  53,  the  combination  of  index  2 
combined  with  index  17  will  provide  the  most  accurate  solution  since  it  has  the 
lowest  condition  number  of  237.9826.  thus,  the  composite  sensitivity  matrix  is 
displayed  in  equation  (4.11)  and  equation  (4.12) 
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To 

test 

the  capability 

of  this  matrix 

,  an 

identical 

FEM 

was  made 

that 

included  damage  in  one  of  the  elements.  Because  we  made  the  sensitivity  matrix 
for  the  flexural  rigidity  (El),  damage  was  simulated  by  decreasing  the  El  value  in 
one  element  by  5%,  and  tested  for  the  length  of  the  beam,  thus  including  all  15 
elements.  As  detailed  in  Chapter  II,  this  is  accomplished  by  solving  equation 
(2.46),  which  has  been  adapted  and  displayed  here  as  equation  (4.13)  and 
equation(4.14): 
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Figures  54  through  Figure  67  show  the  results  of  testing  for  simulated 
damage  in  every  possible  location.  The  simulated  error  is  identified  as  the  blue 
line,  and  solution  provided  using  equation  (4.13)  and  (4.14)  is  in  yellow. 


Figure  54.  Single  Spring  Error  Identification  with  Damage  at  Element  1 
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Figure  55.  Single  Spring  Error  Identification  with  Damage  at  Element  2 


Figure  56.  Single  Spring  Error  Identification  with  Damage  at  Element  3 
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Figure  57.  Single  Spring  Error  Identification  with  Damage  at  Element  4 


Figure  58.  Single  Spring  Error  Identification  with  Damage  at  Element  5 
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Figure  59.  Single  Spring  Error  Identification  with  Damage  at  Element  6 
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Figure  60.  Single  Spring  Error  Identification  with  Damage  at  Element  7 


79 


1000 


-1000 


-2000 


-3000 


-4000 


-5000 


Figure  61 .  Single  Spring  Error  Identification  with  Damage  at  Element  8 


Figure  62.  Single  Spring  Error  Identification  with  Damage  at  Element  9 
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Figure  63.  Single  Spring  Error  Identification  with  Damage  at  Element  10 
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Figure  64.  Single  Spring  Error  Identification  with  Damage  at  Element  1 1 
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Figure  65.  Single  Spring  Error  Identification  with  Damage  at  Element  12 


Figure  66.  Single  Spring  Error  Identification  with  Damage  at  Element  13 
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Figure  67.  Single  Spring  Error  Identification  with  Damage  at  Element  14 


What  Figures  54  through  Figure  67  tell  us  is  the  use  of  the  sensitivity 
matrix  identified  in  equation  (4.13)  will  accurately  find  error  in  all  places  with 
minimal  false  readings. 

E.  SINGLE  MASS  SYSTEM  DYNAMIC  RESPONSE 

Similar  to  the  development  of  the  single  spring  system,  a  single  lumped 
mass  can  be  synthesized  into  the  initial  FRF  of  the  beam  to  provide  new  and 
relevant  information  for  damage  detection  and/or  model  updating.  As  discussed 
before,  the  concept  of  the  idealized,  lumped  mass  allows  us  to  apply  a  numeric 
value  directly  to  the  translational  DOF  of  the  system  without  having  to  add  a 
corresponding  inertial  force. 

1.  Single  Mass  FRF  Synthesis 

Following  the  same  procedure  as  with  the  single  spring  and  the  process 
detailed  in  Chapter  II,  a  lumped  mass  was  synthesized  into  the  FRF  and  checked 
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against  a  physical  modification  of  the  same  structure.  This  check  is  plotted  in 
Figure  68. 


Modified  CSET  FRF  response  check 


|H25  25(n)|f°r  m=0.01294lbm  mass 


Figure  68.  0.227  kg  (0.5  lb)  Lumped  Mass  on  Node  13  FRF  Synthesis 

Check,  H 25,25 


Since  the  FRF  synthesis  and  modified  model  check  lay  directly  on  top  of 
each  other  as  displayed  in  Figure  68,  we  can  infer  that  synthesis  has  been 
correctly  calculated.  Figure  69  shows  the  results  for  a  progressively  increasing 
amount  of  mass. 
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Figure  69.  Lumped  Mass  on  Node  13  FRF  Synthesis  Progression,  H25,25 


As  identified  earlier,  Figure  69  gives  a  detailed  progression  of  the  FRF 
with  the  addition  of  mass.  Unlike  the  synthesis  of  a  spring,  the  beam  maintains 
two  rigid  body  modes.  Additionally,  for  the  point  of  which  the  mass  is  applied,  the 
resonance  peaks  move  down  in  frequency,  with  a  limit  of  the  anti-resonance 
point  from  the  initial  FRF.  Figure  70  displays  the  effects  of  the  added  mass  on 
the  resonance  frequency  for  the  first  nine  modes. 
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Figure  70.  Modal  Frequency  Variance  Due  to  Mass  Addition  at  Node  13 

Figure  70  displays  the  amount  of  change  the  resonance  frequencies 
experience  for  each  mode  as  the  lumped  mass  increases  which  correlates  to  the 
trend  displayed  in  Figure  69.  From  this  graph,  it  appears  there  is  minimal  change 
in  frequency  after  1 .81  kgs  (4  lbs)  of  added  mass. 

2.  Single  Mass  Sensitivity  Analysis 

Following  the  same  process  for  a  single  spring,  to  find  the  values  of  mass 
that  will  provide  that  will  minimize  the  condition  number,  the  values  of  the  highest 
mode  of  the  of  the  sensitivity  matrix  where  graphed  as  the  amount  of  mass 
changed.  Equation  (4.15)  graphically  displays  which  row  is  used,  and  Figures  71 
through  Figure  85  display  the  corresponding  results. 
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Figure  71 .  Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 

Node  2 
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Figure  72.  Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 

Node  3 
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Figure  73.  Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 

Node  4 
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Figure  74.  Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 

Node  5 


Figure  75.  Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 

Node  6 
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Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 
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Figure  77.  Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 

Node  8 
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Figure  78.  Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 

Node  9 


Figure  79.  Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 

Node  10 
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Figure  80.  Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 

Node  1 1 


Figure  81 .  Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 

Node  12 
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Figure  82.  Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 

Node  13 
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Figure  83.  Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 

Node  14 
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Figure  84.  Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 

Node  15 
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Figure  85.  Effects  on  Sensitivity  Value  for  Mass  Addition  Variance  on 

Node  16 
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Figure  71  through  Figure  85  show  the  effects  of  adding  a  lumped  mass  at 
their  respective  nodes.  Similar  to  the  single  spring,  individual  masses  were 
identified  to  maximize  the  condition  number  and  tabulated  in  Table  1 1 . 


Table  1 1 .  Useable  Single  Mass  Additions 


Node 

Added  Mass  kg(lbm) 

2 

0.113(0.25) 

1.814(4.0) 

3 

0.113(0.25) 

1.814(4.0) 

4 

0.113(0.25) 

1.814(4.0) 

5 

0.113(0.25) 

1.814(4.0) 

6 

0.159(0.35) 

1.814(4.0) 

7 

0.113(0.25) 

1.814(4.0) 

8 

0.113(0.25) 

1.814(4.0) 

9 

Not  useable 

10 

0.113(0.25) 

1.814(4.0) 

11 

0.113(0.25) 

1.814(4.0) 

12 

0.159(0.35) 

1.814(4.0) 

13 

0.113(0.25) 

1.814(4.0) 

14 

0.113(0.25) 

1.814(4.0) 

15 

0.113(0.25) 

1.814(4.0) 

16 

0.113(0.25) 

1.814(4.0) 

To  evaluate  the  best  combination  of  sensitivity  matrices,  the  useable 
single  mass  additions  where  placed  into  an  index  as  tabulated  in  Table  12. 
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Table  12.  Single  Mass  Addition  Index  for  Evaluating  Condition  Number 


Index 

Number 

Node 

Mass 

kg(lbm) 

1 

2 

0.113(0.25) 

2 

2 

1.814(4.0) 

3 

3 

0.113(0.25) 

4 

3 

1.814(4.0) 

5 

4 

0.113(0.25) 

6 

4 

1.814(4.0) 

7 

5 

0.113(0.25) 

8 

5 

1.814(4.0) 

9 

6 

0.159(0.35) 

10 

6 

1.814(4.0) 

11 

7 

0.113(0.25) 

12 

7 

1.814(4.0) 

13 

8 

0.113(0.25) 

14 

8 

1.814(4.0) 

Index 

Number 

Node 

Mass 

kg(lbm) 

15 

10 

0.113(0.25) 

16 

10 

1.814(4.0) 

17 

11 

0.113(0.25) 

18 

11 

1.814(4.0) 

19 

12 

0.159(0.35) 

20 

12 

1.814(4.0) 

21 

13 

0.113(0.25) 

22 

13 

1.814(4.0) 

23 

14 

0.113(0.25) 

24 

14 

1.814(4.0) 

25 

15 

0.113(0.25) 

26 

15 

1.814(4.0) 

27 

16 

0.113(0.25) 

28 

16 

1.814(4.0) 

Following  the  same  procedure  as  the  single  spring,  the  two  highest 
baseline  modes  are  combined  with  six  highest  modes  from  two  different  indexes 
and  evaluated  for  the  overall  sensitivity  matrix  condition  number.  The  result  is 
another  matrix  of  784  different  combinations.  Using  the  same  method  as 
described  for  Figure  52,  a  graphic  was  developed  for  the  different  matrix 
condition  numbers  and  displayed  in  Figure  86. 
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Figure  86.  Sensitivity  Matrix  Combination  Condition  Number  Indicator, 

Single  Mass 


In  this  instance,  the  lowest  sensitivity  matrix  condition  number  displayed  in 
Figure  86  occurs  at  index  9,24  with  a  value  of  354.71.  Using  Table  12,  the 
composite  sensitivity  matrix  for  a  single  lumped  mass  should  be: 
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30.602 

40.523 

33.646 

34.824 

34.873 

33.613 

31.827 

30.562 

30.562 

31.827 

33.613 

34.873 

34.824 

33.646 

40.523 

2.373 

5.797 

2.068 

4.947 

3.356 

4.978 

4.387 

4.205 

5.237 

3.372 

5.987 

2.628 

7.574 

3.094 

8.671 

12.506 

10.805 

6.018 

10.048 

3.816 

5.427 

6.783 

2.986 

8.060 

3.307 

6.144 

6.833 

4.751 

9.186 

8.249 

11.494 

9.376 

8.778 

12.799 

11.057 

6.674 

14.046 

9.923 

7.282 

14.310 

9.789 

10.919 

11.537 

7.549 

10.393 

14.866 

16.520 

13.728 

20.894 

22.518 

15.264 

15.292 

22.553 

20.679 

14.785 

22.568 

39.780 

28.048 

28.203 

30.772 

25.515 

21.130 

18.606 

18.458 

20.972 

22.820 

21.557 

18.777 

18.641 

26.454 

48.128 

38.595 

37.976 

43.247 

35.749 

25.291 

25.715 

25.359 

25.411 

25.764 

26.163 

26.262 

26.204 

32.735 

3.064 

7.684 

2.531 

6.171 

3.162 

5.541 

3.873 

4.790 

4.652 

4.053 

5.079 

2.842 

3.922 

1.232 

6.038 

8.298 

7.934 

3.930 

10.022 

4.118 

7.739 

7.689 

4.150 

10.004 

3.957 

8.417 

10.021 

7.544 

6.595 

5.984 

8.680 

4.335 

6.151 

8.455 

3.949 

6.903 

7.977 

3.726 

7.011 

16.233 

26.028 

30.231 

15.351 

10.625 

16.313 

14.528 

9.038 

14.043 

16.703 

10.195 

11.205 

17.280 

11.530 

14.781 

17.061 

17.943 

28.987 

19.365 

21.442 

26.217 

25.068 

20.014 

19.396 

24.236 

26.526 

22.601 

18.144 

21.547 

12.878 

8.853 

43.303 

34.479 

34.430 

34.350 

33.999 

33.712 

33.703 

33.979 

34.324 

34.524 

34.176 

30.995 

14.119 

7.581 

To  test  the  capability  of  this  matrix,  an  identical 


FEM  was  made  that 


included  damage  in  one  of  the  elements.  Because  we  made  the  sensitivity  matrix 


for  the  flexural  rigidity  (El),  damage  was  simulated  by  decreasing  the  El  value  in 
one  element  by  5  percent,  and  tested  for  the  length  of  the  beam,  thus  including 
all  15  elements.  As  detailed  in  Chapter  II,  this  is  accomplished  by  solving 
equation  (2.46),  which  has  been  adapted  and  displayed  here  as  equation  (4.18) 


and  equation  (4.19): 
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(4.18) 


(4.19) 


Figure  87  through  Figure  100  show  the  results  for  testing  for  damage  in 
every  possible  location.  As  with  the  testing  for  a  single  spring  synthesis,  the  blue 
line  indicates  the  actual  damage  of  a  five  percent  reduction  in  El,  and  the  yellow 
is  the  indicated  damage  from  solving  equation  (4.19). 
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Figure  87.  Single  Mass  Error  Identification  with  Damage  at  Element  1 


Figure  88.  Single  Mass  Error  Identification  with  Damage  at  Element  2 
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Figure  89.  Single  Mass  Error  Identification  with  Damage  at  Element  3 


Figure  90.  Single  Mass  Error  Identification  with  Damage  at  Element  4 
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Figure  91 .  Single  Mass  Error  Identification  with  Damage  at  Element  5 


Figure  92.  Single  Mass  Error  Identification  with  Damage  at  Element  6 
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Figure  93.  Single  Mass  Error  Identification  with  Damage  at  Element  7 


Figure  94.  Single  Mass  Error  Identification  with  Damage  at  Element  8 
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Figure  95.  Single  Mass  Error  Identification  with  Damage  at  Element  9 


Figure  96.  Single  Mass  Error  Identification  with  Damage  at  Element  10 
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Element  # 

Figure  97.  Single  Mass  Error  Identification  with  Damage  at  Element  1 1 


Figure  98.  Single  Mass  Error  Identification  with  Damage  at  Element  12 
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Figure  99.  Single  Mass  Error  Identification  with  Damage  at  Element  13 


Figure  100.  Single  Mass  Error  Identification  with  Damage  at  Element  14 
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Figure  87  through  Figure  100  indicate  the  use  of  the  sensitivity  matrix 
identified  in  equation  (4.18)  will  accurately  find  error  in  all  places  with  minimal 
false  readings. 


F.  TWO  SPRING  SYSTEM 

In  previous  thesis  work,  as  described  in  [11],  a  two  pin  solution  was 
developed  concurrently  with  a  single  pin  solution.  For  example,  a  pin  was 
synthesized  onto  nodes  4  and  10  to  provide  new  sensitivity  matrices  and  new 
resonance  frequencies.  In  some  instances,  this  setup  would  provide  better 
damage  detection  capability  than  a  single  pin  solution.  As  such,  this  thesis  also 
explored  this  capability  using  a  two  spring  solution  as  a  proof  of  concept. 

Because  this  thesis  is  exploring  this  concept  as  a  proof  of  concept,  springs 
being  synthesized  onto  the  beam  are  of  equal  value.  Additionally,  this  thesis  will 
continue  to  use  the  same  process  described  for  single  spring  with  a  178.579  kg/ 
mm  (10,000  Ib/in)  spring  as  the  upper  bound  of  the  spring. 

1.  Two  Spring  FRF  Analysis 

The  process  for  synthesizing  a  second  spring  onto  the  beam  follows  the 
same  procedure  described  in  Chapter  II  and  demonstrated  earlier  in  this  chapter. 
The  primary  difference  is  expanding  the  CSET  to  include  the  two  translational 
DOFs,  and  expanding  equation  (2.65)  into  a  2x2  matrix  with  the  spring  values  on 
the  diagonal  so  it  can  be  used  in  equation  (2.64).  As  such,  equation  (2.65)  is 
displayed  here  as  equation  (4.20): 


[AZ(®)] 


AKpring  0 

0  Ak . 


(4.20) 


Figure  101  checks  the  synthesis  of  two  springs  against  a  FEM  that  would 
have  the  springs  physically  installed. 
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Figure  101 .  Two  Springs  at  Nodes  5  and  1 2  FRF  Synthesis  Check,  H33j33, 

178.579  kg/mm  (10,000  Ib/in)  Springs 


Following  the  procedure  in  Chapter  II  with  the  modification  presented  in 
equation  (4.15),  we  can  synthesize  two  springs  onto  the  beam  and  check  the 
synthesis  by  creating  a  second  beam  with  the  new  stiffness  installed.  Using  two 
springs  at  178.579  kg/mm  (10,000  Ib/in)  on  nodes  5  and  12,  Figure  101  confirms 
that  the  springs  are  being  synthesized  correctly. 

As  displayed  in  Figure  102,  when  we  vary  the  strength  of  the  spring  on 
these  nodes,  similar  to  the  single  spring,  the  resonance  frequencies  change, 
confirming  that  useful  information  can  be  used  from  this  process.  Additionally,  we 
see  that  the  addition  of  the  springs  transfers  both  rigid  body  modes  into  elastic 
modes 


107 


500 


450 


400 


—  350 

-S= 

S'  300 
c 
<1) 

o"  250 
a) 

i— 

LL 

ro  200 

"O 

o 

2  150 


100 


50 


0 


0  0.017858  0.17858  0.8929  1.7858  8.929  17.858  178.58  Pin 

Spring  Strength  (kg/mm) 


■Mode  1 
■  Mode  6 


Mode  2 
Mode  7 


-Mode  3 
■Mode  8 


Mode  4 
Mode  9 


Mode  5 


Figure  1 02.  Frequency  Progression  with  Springs  Installed  at  Nodes  5  and 
12,  Varying  strength  from  0-178.579  kg/mm  (0-10,000  Ib/in) 

Equally 


Using  equation  (2.64)  and  (2.65),  we  are  able  to  predict  the  modal 
response  of  the  beam  with  the  spring  installed  as  detailed  in  Figure  103  and 
Figure  104. 
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Figure  103.  FRF  Progression  with  Springs  Installed  at  Nodes  5  and  12, 
Varying  strength  from  0-178.579  kg/mm  (0-10,000  Ib/in) 

Equally,  H99 
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Figure  1 04.  Magnified  FRF  Variance  with  Springs  Installed  at  Nodes  5  and 
12,  Varying  strength  from  0-178.579  kg/mm  (0-10,000  Ib/in) 

Equally,  H9i9 
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2.  Two  Spring  Sensitivity  Analysis 

Following  the  procedure  described  in  the  single  spring  synthesis  analysis, 
this  thesis  explored  the  variance  of  the  highest  mode  of  the  sensitivity  matrix  as 
spring  strength  increased.  A  two-spring  synthesis  can  be  installed  in  105  different 
locations;  however,  as  a  proof  of  concept,  this  thesis  will  not  be  exploring  every 
possible  combination  of  sensitivity  matrices.  As  such,  14  different  locations 
where  chosen  at  random.  These  locations  and  the  subsequent  results  are 
displayed  in  Figure  105  through  Figure  118. 


Figure  1 05.  Effects  on  Sensitivity  Value  with  Increasing  Spring  Strength  for 
Springs  Located  on  Nodes  2  and  6  from  0-178.579  kg/mm 
(0-10,000  Ib/in)  Equally 
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Figure  1 06.  Effects  on  Sensitivity  Value  with  Increasing  Spring  Strength  for 
Springs  Located  on  Nodes  2  and  16  from  0-178.579  kg/mm  (0- 

10,000  Ib/in)  Equally 


Figure  1 07.  Effects  on  Sensitivity  Value  with  Increasing  Spring  Strength  for 
Springs  Located  on  Nodes  3  and  15  from  0-178.579  kg/mm  (0- 

10,000  Ib/in)  Equally 
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Figure  1 08.  Effects  on  Sensitivity  Value  with  Increasing  Spring  Strength  for 
Springs  Located  on  Nodes  4  and  8  from  0-178.579  kg/mm 
(0-10,000  Ib/in)  Equally 


Figure  1 09.  Effects  on  Sensitivity  Value  with  Increasing  Spring  Strength  for 
Springs  Located  on  Nodes  4  and  15  from  0-178.579  kg/mm 
(0-10,000  Ib/in)  Equally 
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Figure  110.  Effects  on  Sensitivity  Value  with  Increasing  Spring  Strength  for 
Springs  Located  on  Nodes  5  and  8  from  0-178.579  kg/mm 
(0-10,000  Ib/in)  Equally 


Figure  111.  Effects  on  Sensitivity  Value  with  Increasing  Spring  Strength  for 
Springs  Located  on  Nodes  5  and  12  from  0-178.579  kg/mm 
(0-10,000  Ib/in)  Equally 
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Figure  112.  Effects  on  Sensitivity  Value  with  Increasing  Spring  Strength  for 
Springs  Located  on  Nodes  5  and  16  from  0-178.579  kg/mm 
(0-10,000  Ib/in)  Equally 


Figure  113.  Effects  on  Sensitivity  Value  with  Increasing  Spring  Strength  for 
Springs  Located  on  Nodes  7  and  1 1  from  0-178.579  kg/mm 
(0-10,000  Ib/in)  Equally 
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gure  114.  Effects  on  Sensitivity  Value  with  Increasing  Spring  Strength  for 
Springs  Located  on  Nodes  7  and  15  from  0-178.579  kg/mm 
(0-10,000  ib/in)  Equally 
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Figure  115.  Effects  on  Sensitivity  Value  with  Increasing  Spring  Strength  for 
Springs  Located  on  Nodes  9  and  13  from  0-178.579  kg/mm 
(0-10,000  Ib/in)  Equally 
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Figure  116.  Effects  on  Sensitivity  Value  with  Increasing  Spring  Strength  for 
Springs  Located  on  Nodes  9  and  16  from  0-178.579  kg/mm 
(0-10,000  Ib/in)  Equally 


Figure  117.  Effects  on  Sensitivity  Value  with  Increasing  Spring  Strength  for 
Springs  Located  on  Nodes  12  and  15  from  0-178.579  kg/mm 
(0-10,000  Ib/in)  Equally 


116 


Figure  118.  Effects  on  Sensitivity  Value  with  Increasing  Spring  Strength  for 
Springs  Located  on  Nodes  15  and  16  from  0-178.579  kg/mm 
(0-10,000  Ib/in)  Equally 

Using  the  14  different  spring  locations  chosen  at  random  as  displayed  in 
Figure  105  through  Figure  1 18  we  can  generate  multiple  sensitivity  matrices  from 
each  set  of  locations.  These  are  identified  in  Table  13. 
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Table  13.  Useable  Two  Spring  Sensitivity  Matrices 


Nodes 

Spring  Value  (kg/mm) 

Spring  Value  (Ib/in) 

2,6 

71.43,  125.01,  178.58 

4000,  7000,  10000 

2,16 

116.08,  178.58 

6500,  10000 

3,15 

53.57,  125.01,  178.58 

3000,  7000,  10000 

4,8 

53.57,  125.01,  178.58 

3000,  7000,  10000 

4,15 

71.43,  142.86,  178.58 

4000,  8000,  10000 

5,8 

53.57,  178.58 

3000,  10000 

5,12 

53.57,  107.15,  125.01,  178.58 

3000,  6000,  7000,  10000 

5,16 

80.36,  178.58 

4500,  1000 

7,11 

53.57,  89.29,  178.58 

3000,  5000,  10000 

7,15 

71.43,  178.58 

4000,  10000 

9,13 

116.08,  160.72 

6500,  9000 

9,16 

125.01,  178.58 

7000,  10000 

12,15 

71.43,  107.15,  160.72,  178.58 

4000,  6000,  9000,  10000 

15,16 

80.36,  178.58 

4500,  1000 

To  minimize  the  sensitivity  matrix  condition  number,  thus  maximizing 
matrix  accuracy,  the  conditions  identified  in  Table  13  are  indexed  and  placed  in 
Table  14.  Following  the  procedures  detailed  for  Figure  52  two  baseline  modes 
are  concatenated  with  six  of  the  highest  modes  from  two  different  indexes.  The 
resulting  sensitivity  matrix  condition  numbers  are  displayed  in  Figure  119. 
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Table  1 4.  Index  of  Two  Spring  Sensitivity  Matrices  for  Evaluating 

Condition  Number 


Index 

Number 

Nodes 

Spring  Value 
kg/mm  (Ib/in) 

1 

2,6 

71.43(4000) 

2 

2,6 

125.01(7000) 

3 

2,6 

178.58(10000) 

4 

2,16 

116.08(6500) 

5 

2,16 

178.58(10000) 

6 

3,15 

53.57(3000) 

7 

3,15 

125.01(7000) 

8 

3,15 

178.58(10000) 

9 

4,8 

53.57(3000) 

10 

4,8 

125.01(7000) 

11 

4,8 

178.58(10000) 

12 

4,15 

71.43(4000) 

13 

4,15 

142.86(8000) 

14 

4,15 

178.58(10000) 

15 

5,8 

53.57(3000) 

16 

5,8 

178.58(10000) 

17 

5,12 

53.57(3000) 

18 

5,12 

107.15(6000) 

Index 

Number 

Nodes 

Spring  Value 
kg/mm  (Ib/in) 

19 

5,12 

125.01(7000) 

20 

5,12 

178.58(10000) 

21 

5,16 

80.36(4500) 

22 

5,16 

178.58(10000) 

23 

7,11 

53.57(3000) 

24 

7,11 

89.29(5000) 

25 

7,11 

178.58(10000) 

26 

7,15 

71.43(4000) 

27 

7,15 

178.58(10000) 

28 

9,13 

116.08(6500) 

29 

9,13 

160.72(9000) 

30 

9,16 

125.01(7000) 

31 

9,16 

178.58(10000) 

32 

12,15 

71.43(4000) 

33 

12,15 

107.15(6000) 

34 

12,15 

160.72(9000) 

35 

12,15 

178.58(10000) 

36 

15,16 

80.36(4500) 

37 

15,16 

178.58(10000) 
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Figure  1 1 9.  Sensitivity  Matrix  Condition  Number  Evaluation,  Two  Springs 


With  37  identified  possibilities,  1,369  sensitivity  matrix  combinations  where 
evaluated  for  their  matrix  condition  number  in  Figure  119.  The  lowest  condition 
number  occurred  at  index  22,30  with  a  condition  number  of  1 12.55.  Using  Table 
14,  the  equation  (4.16)  should  provide  the  most  accurate  combination  of 
conditions  for  solving  for  damage  detection: 
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To  test  the  capability  of  this  matrix,  an  identical  FEM  was  made  that 
included  damage  in  one  of  the  elements.  Because  we  made  the  sensitivity  matrix 
for  the  flexural  rigidity  (El),  damage  was  simulated  by  decreasing  the  El  value  in 


one  element  by  5%,  and  evaluated  all  elements  of  the  beam,  i.e.,  including  all  15 
elements.  As  detailed  in  Chapter  II,  this  is  accomplished  by  solving  equation 


(2.46),  which  has  been  adapted  and  displayed  here  as  equation  (4.23)  and 
equation  (4.24): 
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Figure  120  through  Figure  133  show  the  results  of  testing  for  damage  in 
every  possible  location.  As  with  the  testing  for  a  single  spring  synthesis,  the  blue 
line  indicates  the  actual  damage  of  a  five  percent  reduction  in  El,  and  the  yellow 
is  the  indicated  damage  from  solving  equation  (4.24). 
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Figure  120.  Two  Spring  Error  Identification  at  Element  1 


Figure  1 21 .  Two  Spring  Error  Identification  at  Element  2 
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Figure  122.  Two  Spring  Error  Identification  at  Element  3 


Figure  123.  Two  Spring  Error  Identification  at  Element  4 
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Figure  124.  Two  Spring  Error  Identification  at  Element  5 


Figure  125.  Two  Spring  Error  Identification  at  Element  6 
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Figure  126.  Two  Spring  Error  Identification  at  Element  7 


Figure  127.  Two  Spring  Error  Identification  at  Element  8 
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Figure  128.  Two  Spring  Error  Identification  at  Element  9 


Figure  129.  Two  Spring  Error  Identification  at  Element  10 
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Figure  1 30.  Two  Spring  Error  Identification  at  Element  1 1 


Figure  131.  Two  Spring  Error  Identification  at  Element  12 
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Figure  132.  Two  Spring  Error  Identification  at  Element  13 


Figure  133.  Two  Spring  Error  Identification  at  Element  14 
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Figure  120  through  Figure  133  confirm  that  the  use  of  the  sensitivity  matrix 
identified  in  equation  (4.23)  will  accurately  identify  damage  in  all  locations  with 
minimal  errors 

G.  CONCLUDING  REMARKS 

For  each  scenario  of  synthesis  of  a  single,  dimensionless  spring  to 
ground,  a  single  lumped  mass,  and  two  nodes,  two  springs,  a  combination  of 
conditions  was  identified  that  could  accurately  detect  simulated  damage  in  all 
elements.  The  matrix  condition  number  indicates  the  two  spring  solution  should 
provide  the  most  accurate  solution  when  solving  for  model  updating  and  damage 
detection 

Based  on  the  plots  of  the  sensitivity  value  variance,  there  does  not  appear 
to  be  any  more  information  provided  by  adding  a  single  lumped  mass.  However, 
further  exploration  should  be  conducted  in  the  exploration  of  the  single  spring 
and  two  spring  scenarios  to  discover  their  limitations  and  capability  to  add  more 
information  to  the  sensitivity  matrix. 
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V.  MODEL  UPDATING  AND  DAMAGE  DETECTION 


A.  INTRODUCTION 

As  described  in  Chapter  II,  The  purpose  of  FRF  synthesis  is  to  manipulate 
test  data  to  determine  the  exact  dynamic  response  of  the  beam  due  to  the 
integration  of  a  boundary  condition  or  beam  modification.  The  advantage  of 
which  is  the  ability  to  determine  the  new  resonance  frequencies  without  the 
tedious  implications  of  physically  modifying  the  beam.  These  new  resonance 
frequencies  can  then  be  compared  against  the  correlating  FEM  frequencies  and 
used  in  damage  detection  and  model  updating. 

Experimental  data  will  yield  FRFs,  and  by  using  MATLAB’s  internal 
findpeaks  function  [17],  the  FRFs  will  provide  resonance  frequencies  of  the 
prototype  structure.  Additionally,  by  following  the  procedure  for  FRF  synthesis 
detailed  in  Chapter  II,  we  can  synthesize  the  beam  modifications  demonstrated  in 
Chapter  IV  on  measured  data  to  update  the  FEM  and  perform  damage  detection 

It  is  important  to  note  that  node  numbering  for  the  FEM  are  not  an  exact 
match  for  the  nodes  on  the  measured  beam.  For  example,  node  2  on  the  FEM  is 
physical  node  1  on  the  measured  beam.  Additionally,  the  data  captured  from  the 
measured  beam  captures  the  translational  movement  only.  As  such,  the  DOF 
number  for  the  measured  beam  reflects  the  node  number  directly. 

B.  VERIFICATION  OF  DATA  IMPORT 

All  experimental  data  gathered  in  Chapter  III  using  the  Single  Impact 
Hammer  test  was  exported  to  Microsoft  Excel  .xlsx  files  and  imported  into  Matlab 
using  the  Matlab  proprietary  xlsread  function  [17].  By  plotting  the  driving  points, 
we  can  confirm  the  data  was  successfully  transferred  by  ensuring  the  resonance 
peak  frequencies  line  up  as  displayed  in  Figure  134.  If  the  peaks  did  not  occur  at 
the  same  frequency,  we  would  suspect  an  issue  with  either  the  measurement  or 
data  import. 
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Figure  134.  Experimental  Beam  Response  Driving  Point  Data  Verification 

C.  MODEL  UPDATING  USING  SINGLE  SPRING 

Following  the  procedure  detailed  in  Chapter  II  and  demonstrated  in 
Chapter  IV,  the  conditions  detailed  in  equation  (4.13)  were  synthesized  into  the 
experimental  data.  First,  a  178.579  kg/mm  (10,000  Ib/in)  spring  was  synthesized 
onto  node  1  using  equation  (2.64)  and  provides  the  dynamic  response  detailed  in 
Figure  135. 
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Figure  135.  Synthesized  Measured  Data,  178.579  kg/mm  (10,000  Ib/in) 

Spring  Synthesized  on  Node  1 

In  order  to  solve  equation  (4.13),  we  need  to  determine  the  new 
resonance  frequencies  of  the  beam.  To  do  so,  all  driving  point  FRFs  are 
averaged,  providing  a  single  FRF  as  displayed  in  Figure  136.  Additionally,  this 
minimizes  false  positives  in  identifying  resonance  frequencies  and  provides  a 
clear  comparison  to  the  averaged  FRFs  of  the  FEM. 
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Figure  136.  Synthesized  Measured  Data,  178.579  kg/mm  (10,000  Ib/in)  on 
Node  1,  Averaged  Driving  Points  compared  against  FEM 


Using  the  MATLAB  findpeaks  function  [17],  we  can  determine  the 
resonance  frequencies  of  the  synthesized  measured  data  and,  as  determined  in 
the  development  of  equation  (4.13),  use  the  six  highest  modes  in  damage 
detection. 

To  get  the  remaining  modes,  we  synthesize  a  178.579  kg/mm  (10,000  lb/ 
in)  spring  onto  the  baseline  FRFs  at  node  7  using  equation  (2.64). 
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Figure  137.  Synthesized  Measured  Data,  178.579  kg/mm  (10,000  Ib/in) 

Spring  Synthesized  on  Node  7 

Synthesis  of  a  178.579  kg/mm  (10,000  Ib/in)  spring  onto  the  baseline 
FRFs  at  node  7  produces  0,  when  looking  at  the  driving  points.  Following  the 
same  procedure  as  the  single  spring,  all  driving  point  FRFs  are  averaged, 
providing  a  single  FRF  as  displayed  in  Figure  138. 
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Figure  138.  Synthesized  Measured  Data,  178.579  kg/mm  (10,000  Ib/in)  on 
Node  7,  Averaged  Driving  Points  compared  against  FEM 


Again,  using  the  MATLAB  findpeaks  function  [17],  we  can  determine  the 
resonance  frequencies  of  the  synthesized  measured  data  and,  as  determined  in 
the  development  of  equation  (4.8),  use  the  six  highest  modes  in  damage 

detection.  Understanding  that  |Ao2}  =  {A/l}  where  go  is  the  resonance  frequency, 

the  difference  in  resonance  frequencies  between  the  six  modes  identified  in  the 
measured  data  and  the  synthesized  FEM  are  then  concatenated  to  the  six 
identified  with  the  previous  synthesis  and  two  of  the  baseline  modes  to  create  the 

{A/l}  vector  identified  in  equation  (4.13)  and  reflected  here  as  equation  (5.1), 
modified  to  reflect  the  synthesized  nodes  for  the  measured  data. 
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Using  the  sensitivity  matrix  identified  for  a  single  spring  in  Chapter  IV  and 
the  difference  in  resonance  frequencies  developed  using  synthesis,  the  solution 
to  equation(5.1)  is  displayed  in  Figure  139. 


Figure  139.  Identified  Error  in  FEM  Flexural  Rigidity  Using  Single  Spring 

Synthesis 
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The  results  of  equation  (5.1),  displayed  as  Figure  139,  equations  (5.2)  and 
equations  (5.3),  provide  the  necessary  corrections  to  the  FEM  elements  to  better 
reflect  the  prototype’s  dynamic  response.  Using  these  corrections,  the  FRF  for 
the  FEM  was  recalculated  using  the  element  specific  correction  identified  in 
equation  (5.2),  displayed  in  Figure  140  with  Table  15.  Separately,  the  FEM  was 
recalculated  using  the  average  element  correction  identified  in  equation  (5.3)  and 
displayed  in  Figure  141  with  Table  16. 
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Figure  140.  Updated  FEM  FRF  compared  to  Prototype  FRF  Using  Element 

Specific  Correction,  Flu 


Table  15.  Updated  FEM  Resonance  Frequencies  Using  Single  Spring 

Synthesis  Element  Specific  Correction 


Previous  FEM 
Frequency(Hz) 

Previous 

%  Rel 

Error 

Updated 
FEM  Freq 
(Hz) 

Measured 
Freq  (Hz) 

Absolute 

Error 

%  Relative 

Error 

%  Error 

Corrected 

14.13452 

2.7965 

13.996 

13.75 

0.246 

1.7891 

1.00742 

38.95759 

2.1838 

36.7333 

38.125 

-1.3917 

-3.6504 

-1.46652 

76.37207 

1.8294 

74.5518 

75 

-0.4482 

-0.5976 

1.23183 

126.27156 

1.7806 

124.5348 

124.0625 

0.4723 

0.3807 

1.39991 

188.72298 

2.0124 

186.5883 

185 

1.5883 

0.8585 

1.15388 

263.83435 

2.0883 

261.7421 

258.4375 

3.3046 

1.2787 

0.80958 

351.78075 

2.3362 

348.7656 

343.75 

5.0156 

1.4591 

0.87713 

452.81294 

2.5479 

448.5534 

441.5625 

6.9909 

1.5832 

0.96465 

567.24773 

2.8438 

563.0675 

551.5625 

11.505 

2.0859 

0.75789 

695.40391 

3.4058 

687.8183 

672.5 

15.3183 

2.2778 

1.12797 

837.35461 

3.8982 

828.9641 

805.9375 

23.0266 

2.8571 

1.04109 

991.88069 

4.3056 

980.0004 

950.9375 

29.0629 

3.0562 

1.24932 
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Updated  FEM  Using  Average  Spring  Correction  vs  Measured  Data 
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Figure  141 .  Updated  FEM  FRF  Compared  to  Prototype  FRF  Using  Average 

Element  Correction 


Table  16.  Updated  FEM  Resonance  Frequencies  Using  Single  Spring 

Synthesis  Average  Element  Correction 


Previous  FEM 
Frequency(Hz) 

Previous 

%  Rel 

Error 

Updated 
FEM  Freq 
(Hz) 

Measured 
Freq  (Hz) 

Absolute 

Error 

%  Relative 

Error 

%  Error 

Corrected 

14.13452 

2.7965 

14.1235 

13.75 

0.3735 

2.7164 

0.08015 

38.95759 

2.1838 

38.9271 

38.125 

0.8021 

2.1039 

0.07997 

76.37207 

1.8294 

76.3123 

75 

1.3123 

1.7497 

0.07969 

126.27156 

1.7806 

126.1728 

124.0625 

2.1103 

1.7010 

0.07961 

188.72298 

2.0124 

188.5753 

185 

3.5753 

1.9326 

0.07983 

263.83435 

2.0883 

263.6279 

258.4375 

5.1904 

2.0084 

0.07988 

351.78075 

2.3362 

351.5055 

343.75 

7.7555 

2.2561 

0.08007 

452.81294 

2.5479 

452.4587 

441.5625 

10.8962 

2.4676 

0.08022 

567.24773 

2.8438 

566.8039 

551.5625 

15.2414 

2.7633 

0.08047 

695.40391 

3.4058 

694.8599 

672.5 

22.3599 

3.3249 

0.08089 

837.35461 

3.8982 

836.6995 

805.9375 

30.762 

3.8169 

0.08129 

991.88069 

4.3056 

991.1047 

950.9375 

40.1672 

4.2240 

0.08160 
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D.  MODEL  UPDATING  USING  SINGLE  MASS 

Following  the  same  procedures  described  for  the  synthesis  of  a  single 
spring,  the  conditions  of  the  lumped  mass  identified  in  equation  (4.17)  where 
synthesized  into  the  Measured,  baseline  FRF.  Equation  (4.19)  is  restated  here 
with  the  modification  to  node  locations  to  reflect  placement  on  the  beam. 


S(0) 

J2xl4 

'  ^(N 5-035)' 

J6xl4 

£(#13-4.00)" 

J6xl4 


"\ 

o 

_ J 

A  El, 

t  )  2x1 

|a2("1 * * 4 5 * 7 8 * 10 * *-0  35)1 

>  =  < 

aei2 

> 

l  )  6x1 

I 

Ia/I'"13 14 15-4'0”'! 

A  J6xl, 

A  EIU 

(5.4) 


First,  a  mass  of  0.1588  kg  (0.35  lb)  was  synthesized  onto  node  5.  The 
resulting  synthesized  FRFs  generated  are  displayed  in  Figure  142. 
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Figure  142.  Synthesized  Measured  Data,  0.1588  kg  (0.35  lb)  on  Node  5 
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By  averaging  the  driving  points  identified  in  Figure  142,  we  can  once  again 
identify  the  resonance  frequencies  required  in  both  the  synthesized  measured 
data  and  synthesized  FEM  FRFs,  as  displayed  in  Figure  143,  for  inclusion  in 
equation  (5.4). 
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Figure  143.  Synthesized  Measured  Data,  Averaged  Driving  Points  0.1588 

kg  (0.35  lb)  on  Node  5 

Next,  a  mass  of  1 .81437kg  (41b)  is  synthesized  into  the  measured  data  on 
node  13.  The  resulting  synthesized  FRFs  are  displayed  as  Figure  144  With 
driving  point  averages  displayed  as  Figure  145. 
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Figure  144.  Synthesized  Measured  Data,  1.81437kg  (41b)  on  Node  13 


Figure  145.  Synthesized  Measured  Data,  Average  Driving  Points, 

1 ,81437kg  (41b)  on  Node  13 
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Again,  using  the  MATLAB  findpeaks  function  [17]  on  Figure  145  as 
averaged  from  Figure  144  we  can  determine  the  resonance  frequencies  of  the 
synthesized  measured  data  and  the  synthesized  FEM.  As  determined  in  the 
development  of  equation  (5.4)  the  differences  in  the  six  highest  modes  are  used 
in  damage  detection. 

In  using  the  sensitivity  matrix  identified  in  equation  (4.19)  and  finding  the 
differences  in  resonance  frequencies  of  the  required  modes  as  identified  by 
equation  (5.4),  Figure  143,  and  Figure  145,  we  can  solve  for  the  errors  in  the 
FEM  flexural  rigidity  vector  as  displayed  in  Figure  146. 


Figure  146.  Identified  Error  in  FEM  Flexural  Rigidity  using  Single  Mass 

Synthesis 
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The  results  of  equation  (5.4),  displayed  in  Figure  146,  equation  (5.5),  and 
equation  (5.6),  provide  the  necessary  corrections  to  the  FEM  elements  to  better 
reflect  the  prototype’s  dynamic  response.  Using  these  corrections,  the  FRF  for 
the  FEM  was  recalculated  using  the  element  specific  correction  identified  in 
equation  (5.5),  displayed  in  Figure  147  with  Table  17.  Separately,  the  FEM  was 
recalculated  using  the  average  element  correction  identified  in  equation  (5.6)  and 
displayed  in  Figure  148  with  Table  18. 
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Figure  147.  Updated  FEM  FRF  compared  to  Prototype  FRF  Using  Mass 

Synthesis  Element  Specific  Correction,  Flu 


Table  17.  Updated  FEM  Resonance  Frequencies  Using  Mass  Synthesis 

Element  Specific  Correction 


Previous  FEM 

Previous 

%  Rel 

Error 

Updated  FEM 

Measured 

Absolute 

%  Relative 

%  Error 

Frequency(Hz) 

Freq  (Hz) 

Freq  (Hz) 

Error 

Error 

Corrected 

14.13452 

2.3291 

13.6115 

13.81281 

-0.20127 

-1.4571 

0.87196 

38.95759 

1.7872 

35.7128 

38.27358 

-2.56082 

-6.6908 

-4.90367 

76.37207 

1.7888 

71.4144 

75.02992 

-3.61556 

-4.8188 

-3.03000 

126.27156 

1.8210 

118.0948 

124.01327 

-5.91842 

-4.7724 

-2.95141 

188.72298 

1.8938 

184.9390 

185.21543 

-0.2764 

-0.1492 

1.74454 

263.83435 

2.0310 

257.2172 

258.58254 

-1.36536 

-0.5280 

1.50298 

351.78075 

2.1856 

344.0446 

344.25651 

-0.2119 

-0.0616 

2.12410 

452.81294 

2.3864 

443.6277 

442.25887 

1.368806 

0.3095 

2.07690 

567.24773 

2.6374 

553.5058 

552.67146 

0.834297 

0.1510 

2.48646 

695.40391 

3.0958 

679.6416 

674.52174 

5.119838 

0.7590 

2.33682 

837.35461 

3.5819 

815.3904 

808.39895 

6.991492 

0.8649 

2.71700 

991.88069 

3.9774 

954.7796 

953.93877 

0.840799 

0.0881 

3.88926 
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Figure  148.  Updated  FEM  FRF  compared  to  Prototype  FRF  Using  Mass 

Synthesis  Average  Element  Correction, 


Table  18.  Updated  FEM  Resonance  Frequencies  Using  Mass  Synthesis 

Average  Element  Correction 


Previous  FEM 

Previous 

%  Rel 

Error 

Updated  FEM 

Measured 

Absolute 

%  Relative 

%  Error 

Frequency(Hz) 

Freq  (Hz) 

Freq  (Hz) 

Error 

Error 

Corrected 

14.13452 

2.3291 

14.1124 

13.81281 

0.299551 

2.1686 

0.16042 

38.95759 

1.7872 

38.8965 

38.27358 

0.622948 

1.6276 

0.15954 

76.37207 

1.7888 

76.2524 

75.02992 

1.222447 

1.6293 

0.15954 

126.27156 

1.8210 

126.0736 

124.01327 

2.060375 

1.6614 

0.15959 

188.72298 

1.8938 

188.4272 

185.21543 

3.211751 

1.7341 

0.15971 

263.83435 

2.0310 

263.4208 

258.58254 

4.838285 

1.8711 

0.15992 

351.78075 

2.1856 

351.2294 

344.25651 

6.972875 

2.0255 

0.16016 

452.81294 

2.3864 

452.1032 

442.25887 

9.844347 

2.2259 

0.16048 

567.24773 

2.6374 

566.3587 

552.67146 

13.6872 

2.4766 

0.16087 

695.40391 

3.0958 

694.3140 

674.52174 

19.79222 

2.9343 

0.16159 

837.35461 

3.5819 

836.0422 

808.39895 

27.64324 

3.4195 

0.16235 

991.88069 

3.9774 

990.3261 

953.93877 

36.3873 

3.8144 

0.16297 
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E.  MODEL  UPDATING  USING  TWO  SPRINGS 

Following  the  same  procedures  described  for  the  synthesis  of  a  single 
spring,  two  springs  of  equal  strength,  as  identified  in  equation  (4.21),  were 
synthesized  into  the  measured,  baseline  FRF.  Equation  (4.23)  is  restated  here 
with  the  modification  to  node  locations  to  reflect  placement  on  the  beam. 
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First,  two  springs,  each  of  strength  178.579  kg/mm  (10,000  Ib/in)  were 
synthesized  onto  node  4  and  1 1 : 
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Figure  149.  Synthesized  Measured  Data,  178.579  kg/mm  (10,000  Ib/in)  on 

Nodes  4  and  1 1 
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Figure  150.  Synthesized  Measured  Data,  Averaged  Driving  Points  178.579 

kg/mm  (10,000  Ib/in)  on  Nodes  4  and  1 1 

By  averaging  the  driving  points  identified  in  Figure  149,  we  can  once  again 
identify  the  resonance  frequencies  required  in  both  the  synthesized  measured 
data  and  synthesized  FEM  FRFs,  as  displayed  in  Figure  150,  for  inclusion  in 
equation  (5.7). 

Next  two  springs,  each  of  strength  178.579  kg/mm  (10,000  Ib/in)  were 
synthesized  onto  node  4  and  15: 
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Figure  151.  Synthesized  Measured  Data,  178.579  kg/mm  (10,000  Ib/in)  on 

Nodes  4  and  15 


Figure  152.  Synthesized  Measured  Data,  Average  Driving  Points,  1178.579 

kg/mm  (10,000  Ib/in)  on  Nodes  4  and  15 
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Again,  using  the  MATLAB  findpeaks  function  [17]  on  Figure  152,  as 
averaged  from  Figure  151,  we  can  determine  the  resonance  frequencies  of  the 
synthesized  measured  data  and  the  synthesized  FEM.  As  determined  in  the 
development  of  equation  (5.4)  the  differences  in  the  six  highest  modes  are  used 
in  damage  detection. 

In  using  the  sensitivity  matrix  identified  in  equation  (4.12)  and  finding  the 
differences  in  resonance  frequencies  of  the  required  modes  as  identified  by 
equation  (5.7),  Figure  150,  and  Figure  152,  we  can  solve  for  the  errors  in  the 
FEM  flexural  rigidity  vector: 


1  2  3  4  5  6  7  8  9  10  11  12  13  14 

Element# 


Figure  153.  Identified  Error  in  FEM  Flexural  Rigidity  Using  Two  Spring 

Synthesis 
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{AEI} 


'3413.727 

-32059.664 

18218.355 

-10847.512 

24591.878 

12268.331 

36449.882 

13812.577 

-20482.147 

-39399.042 

-17706.379 

-11632.960 

7529.408 

14729.118 


{AEI} 

AVERAGE 


-79.602 


(5.8) 


(5.9) 


The  results  of  equation  (5.7),  displayed  in  Figure  153,  equation  (5.8),  and 
equation  (5.9),  provide  the  necessary  corrections  to  the  FEM  elements  to  better 
reflect  the  prototype’s  dynamic  response.  Using  these  corrections,  the  FRF  for 
the  FEM  was  recalculated  using  the  element  specific  correction  identified  in 
equation  (5.8),  displayed  in  Figure  154  with  Table  19.  Separately,  the  FEM  was 
recalculated  using  the  average  element  correction  identified  in  equation  (5.9)  and 
displayed  in  Figure  155  with  Table  20. 
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Figure  154.  Updated  FEM  FRF  compared  to  Prototype  FRF  Using  Two 
Spring  Synthesis  Element  Specific  Correction,  Flu 


Table  19.  Updated  FEM  Resonance  Frequencies  Using  Two  Spring 

Synthesis  Element  Specific  Correction 


Previous  FEM 

Previous 

%  Rel 

Error 

Updated 

FEM  Freq 
(Hz) 

Measured 

Absolute 

%  Relative 

%  Error 

Frequency(Hz) 

Freq  (Hz) 

Error 

Error 

Corrected 

14.13452 

2.3291 

13.7422 

13.81281 

-0.07057 

-0.5109 

1.81816 

38.95759 

1.7872 

36.4948 

38.27358 

-1.77878 

-4.6475 

-2.86038 

76.37207 

1.7888 

74.8184 

75.02992 

-0.21157 

-0.2820 

1.50684 

126.27156 

1.8210 

120.4998 

124.01327 

-3.51352 

-2.8332 

-1.01217 

188.72298 

1.8938 

179.6513 

185.21543 

-5.56412 

-3.0041 

-1.11037 

263.83435 

2.0310 

249.6324 

258.58254 

-8.95011 

-3.4612 

-1.43022 

351.78075 

2.1856 

340.3857 

344.25651 

-3.87084 

-1.1244 

1.06124 

452.81294 

2.3864 

443.0534 

442.25887 

0.794519 

0.1797 

2.20675 

567.24773 

2.6374 

554.6007 

552.67146 

1.929261 

0.3491 

2.28834 

695.40391 

3.0958 

679.7198 

674.52174 

5.198101 

0.7706 

2.32521 

837.35461 

3.5819 

815.7007 

808.39895 

7.301723 

0.9032 

2.67862 

991.88069 

3.9774 

960.9866 

953.93877 

7.047855 

0.7388 

3.23858 
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Figure  155.  Updated  FEM  FRF  compared  to  Prototype  FRF  Using  Two 
Spring  Synthesis  Average  Element  Correction,  Hi,i 


Table  20.  Updated  FEM  Resonance  Frequencies  Using  Two  Spring 

Synthesis  Average  Element  Correction 


Previous  FEM 
Frequency(Hz) 

Previous 

%  Rel 

Error 

Updated 

FEM  Freq 
(Hz) 

Measured 
Freq  (Hz) 

Absolute 

Error 

%  Relative 

Error 

%  Error 

Corrected 

14.13452 

2.3291 

14.1276 

13.81281 

0.314755 

2.2787 

0.05035 

38.95759 

1.7872 

38.9384 

38.27358 

0.664852 

1.7371 

0.05005 

76.37207 

1.7888 

76.3345 

75.02992 

1.304595 

1.7388 

0.05005 

126.27156 

1.8210 

126.2095 

124.01327 

2.196196 

1.7709 

0.05007 

188.72298 

1.8938 

188.6302 

185.21543 

3.414747 

1.8437 

0.05011 

263.83435 

2.0310 

263.7046 

258.58254 

5.122073 

1.9808 

0.05017 

351.78075 

2.1856 

351.6078 

344.25651 

7.35126 

2.1354 

0.05025 

452.81294 

2.3864 

452.5903 

442.25887 

10.3314 

2.3361 

0.05035 

567.24773 

2.6374 

566.9688 

552.67146 

14.29734 

2.5870 

0.05047 

695.40391 

3.0958 

695.0620 

674.52174 

20.54022 

3.0452 

0.05070 

837.35461 

3.5819 

836.9429 

808.39895 

28.54391 

3.5309 

0.05093 

991.88069 

3.9774 

991.3930 

953.93877 

37.45419 

3.9263 

0.05113 
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F.  MODEL  UPDATING  USING  LOWER  MODES 

The  results  from  model  updating  with  the  single  spring,  single  mass,  and 
two  spring  solutions  all  utilized  the  highest  six  modes  from  each  index  in  Table 
10,  Table  12,  and  Table  14  respectively.  Interestingly  enough,  the  highest  six 
modes  in  the  base  FEM  after  updating  have  the  best  error  correction.  Because  of 
this,  a  short  test  was  conducted  using  the  two  spring  setup  and  conditions 
identified  in  Table  14,  but  including  some  lower  modes.  As  such,  the  new 
composite  sensitivity  matrix  was  made  as  follows: 


^  (A/4,1 1-10,000) 


2,3,4,5,11,12,13 


\s  1  = 

[_  comp  J 


L 

comp  _| 

[14x14 

r 
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LL 

J6,7,8,9, 10, 12,13 

- 
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0.000 

0.008 

0.037 

0.056 

0.044 

0.027 

0.010 

0.001 

0.007 

0.032 

0.035 

0.013 

0.003 

0.000“ 

0.003 

0.055 

0.235 

0.185 

0.011 

0.067 

0.149 

0.112 

0.019 

0.064 

0.132 

0.060 

0.015 

0.001 

0.014 

0.164 

0.437 

0.139 

0.665 

0.762 

0.095 

0.573 

0.996 

0.247 

0.317 

1.098 

0.723 

0.075 

(5.11) 

0.017 

0.163 

0.296 

0.114 

0.676 

0.393 

0.140 

0.667 

0.210 

1.144 

0.752 

1.831 

2.159 

0.278 

25.431 

18.162 

24.307 

11.328 

7.439 

12.584 

11.584 

6.903 

11.534 

6.879 

2.064 

8.072 

4.845 

6.766 

9.239 

6.057 

8.434 

13.133 

10.829 

9.297 

12.989 

14.299 

11.124 

8.846 

28.775 

32.559 

24.815 

38.110 

19.713 

14.569 

13.363 

25.734 

33.431 

30.279 

27.885 

28.695 

31.873 

30.058 

21.234 

18.141 

18.294 

25.347 

0.116 

0.558 

0.145 

0.340 

0.256 

0.334 

1.395 

1.539 

4.644 

0.985 

5.708 

1.460 

4.307 

2.901 

2.168 

8.441 

1.553 

6.928 

2.102 

5.066 

2.155 

1.217 

0.118 

0.060 

0.245 

0.130 

0.179 

0.159 

0.218 

0.419 

0.236 

0.262 

0.390 

0.248 

0.369 

6.339 

4.251 

12.848 

3.801 

12.856 

4.470 

8.920 

9.714 

14.249 

12.462 

6.978 

16.758 

5.473 

10.312 

1.466 

1.625 

1.138 

0.721 

1.612 

0.530 

1.144 

4.408 

3.555 

5.485 

3.375 

3.158 

6.100 

0.688 

12.001 

12.386 

11.240 

19.563 

12.537 

10.988 

12.415 

13.403 

9.041 

9.765 

11.857 

11.652 

9.933 

7.836 

19.783 

26.144 

21.081 

19.607 

23.600 

29.046 

4.069 

39.773 

30.647 

30.104 

30.654 

31.396 

32.259 

23.952 

14.793 

17.195 

17.217 

17.606 

17.455 

19.028 

6.485 

This  composite  sensitivity  matrix  does  not  include  any  baseline  modes; 
however,  it  includes  modes  2,  3,  4,  5,  11,  12, and  13  from  index  20  and  modes 
6,7,8,9,10,12,  and  13  from  index  31  indicated  in  Table  14.  The  condition  number 
for  this  matrix  is  4781.7.  This  matrix  conditioning  is  not  great,  but  within  an 
acceptable  magnitude  demonstrated  during  previous  tests. 


Figure  156  and  Figure  157  demonstrate  the  resulting  synthesized  FRFs  of 
the  measured  data  for  a  178.579  kg/mm  (10,000  Ib/in)  spring  synthesized  on 
Nodes  4  and  11.  Similarly,  Figure  158  and  Figure  159  demonstrate  the  resulting 
synthesized  FRFs  of  the  measured  data  for  a  178.579  kg/mm  (10,000  Ib/in) 
spring  synthesized  on  Nodes  8  and  15. 
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Figure  156.  Synthesized  Measured  Data,  178.579  kg/mm  (10,000  Ib/in)  on 

Nodes  4  and  1 1 
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Figure  157.  Synthesized  Measured  Data,  Driving  Point  Average,  178.579 

kg/mm  (10,000  Ib/in)  on  Nodes  4  and  1 1 
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Figure  158.  Synthesized  Measured  Data,  178.579  kg/mm  (10,000  Ib/in)  on 

Nodes  8  and  15 


Figure  159.  Synthesized  Measured  Data,  Driving  Point  Average,  178.579 

kg/mm  (10,000  Ib/in)  on  Nodes  8  and  15 
157 


Figure  156  through  Figure  159  demonstrate  the  synthesis  of  the  two 
spring  system  for  the  conditions  indicated  in  equation  (5.10).  Once  synthesized, 
we  extract  the  frequencies  for  the  specific  modes  for  model  updating.  By 
modifying  the  sensitivity  matrix  as  such,  equation  (4.23)  can  be  solved  with  the 
following  modifications: 
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(5.12) 


The  identified  modifications  to  the  El  vector  by  solving  equation  (5.12)  are 
shown  in  Figure  160. 


Figure  160.  Identified  Error  in  FEM  Flexural  Rigidity  Using  Two  Spring 

Synthesis,  Low  Modes 


158 


{A£7}  =  < 
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-1422.751 

13023.016 

4826.005 
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-2920.809 

-17305.534 

13729.188 

-323.949 

-711.184 

-13392.979 

-16949.322 

12155.394 

18278.053 


(5.13) 


lA^L™=-767-69  (5-14) 

The  modification  to  the  El  vector  as  displayed  in  Figure  160  was 
calculated  as  equation  (5.13)  and  averaged  in  equation  (5.14).  By  using  the 
element  specific  correction,  we  can  modify  the  dynamic  response  of  the  beam  as 
calculated  for  Figure  161  and  Table  21 . 
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Figure  161.  Updated  FEM  FRF  compared  to  Prototype  FRF  Using  Two 
Spring  Synthesis  Element  Specific  Correction,  Low  Modes,  Hi,i 


Table  21 .  Updated  FEM  Resonance  Frequencies  Using  Low  Modes 

Synthesis  Element  Specific  Correction 


Previous  FEM 
Frequency(Hz) 

Previous 

%  Rel 

Error 

Updated 
FEM  Freq 
(Hz) 

Measured 
Freq  (Hz) 

Absolute 

Error 

% 

Relative 

Error 

%  Error 

Corrected 

14.1345 

2.3290 

13.9022 

13.8128 

-0.0894 

-0.6474 

1.6816 

38.9576 

1.7872 

38.3917 

38.2736 

-0.1181 

-0.3086 

1.4786 

76.3721 

1.7888 

75.1602 

75.0299 

-0.1302 

-0.1736 

1.6152 

126.2716 

1.8210 

125.4511 

124.0133 

-1.4379 

-1.1594 

0.6616 

188.7230 

1.8938 

188.6115 

185.2154 

-3.3960 

-1.8336 

0.0602 

263.8343 

2.0310 

263.3370 

258.5825 

-4.7545 

-1.8387 

0.1923 

351.7808 

2.1856 

351.5785 

344.2565 

-7.3219 

-2.1269 

0.0588 

452.8129 

2.3864 

454.3128 

442.2589 

-12.0539 

-2.7255 

-0.3391 

567.2477 

2.6374 

566.7458 

552.6715 

-14.0743 

-2.5466 

0.0908 

695.4039 

3.0958 

694.1887 

674.5217 

-19.6669 

-2.9157 

0.1802 

837.3546 

3.5819 

836.4128 

808.3990 

-28.0139 

-3.4654 

0.1165 

991.8807 

3.9774 

991.8157 

953.9388 

-37.8770 

-3.9706 

0.0068 
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G.  CONCLUDING  REMARKS 

As  had  been  developed  in  Chapter  IV,  we  were  able  to  synthesize  a  single 
spring,  single  mass,  and  two  spring  system  into  the  measured  data  and 
successfully  update  the  original  FEM.  Interestingly,  the  high  modes  where 
corrected  almost  exactly  when  using  the  element  specific  corrections;  however, 
low  modes  did  not  see  the  same  amount  of  correction,  or  got  worse.  This  may 
have  been  due  to  the  use  of  high  modes  only  when  developing  the  sensitivity  for 
each  synthesis  type.  The  correction  to  the  low  modes  did  occur  when  the  low 
modes  where  included  into  the  sensitivity  matrix,  but  did  not  provide  a  significant 
correction  to  the  high  modes. 

Overall,  the  Single  Mass  synthesis  provided  the  greatest  amount  of 
correction,  especially  at  the  higher  modes,  but  had  the  largest  increase  in  error  in 
the  low  modes.  The  Two  Spring  synthesis  provided  the  greatest  benefit  with  least 
amount  of  additional  error.  Finally,  the  Two  Spring  Synthesis  with  low  modes 
provided  the  least  amount  of  added  error  but  was  minor  in  the  amount  of 
correction. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

As  identified  in  Chapter  I,  the  primary  purpose  of  this  thesis  was  to  determine 
the  possibility  and  relevance  of  using  an  artificial  spring  to  ground  or  lumped  mass 
added  to  a  beam  in-lieu  of  a  pin  connection.  In  three  different  types  of  simulations, 
we  gathered  new  and  relevant  information  for  damage  detection  and  model  updating 
by  varying  the  magnitude  of  either  a  single  spring,  single  mass,  or  two  spring 
system. 

As  was  theorized,  by  increasing  the  stiffness  of  an  artificial  spring  to  ground, 
the  beam’s  dynamic  response  will  change  to  reflect  that  of  a  pin  connection; 
however,  there  is  a  lot  of  useable  FRF  and  sensitivity  information  in-between  the 
baseline  beam  with  no  spring,  and  a  beam  with  a  spring  of  infinite  stiffness.  By 
varying  the  magnitude  of  the  springs  or  mass  artificially  attached,  we  gathered  new 
frequencies  and  sensitivity  matrices,  which  we  used  in  creating  a  new,  composite 
sensitivity  matrix. 

Following  previous  thesis  work  [11],  high  modes  from  the  new  sensitivity 
matrices  where  developed  into  a  square,  full  rank,  composite  matrix.  Flowever,  due 
to  the  nature  of  the  sensitivity  matrix,  it  was  discovered  the  matrix  condition  number 
would  dramatically  affect  the  results  of  any  damage  detection.  As  such,  multiple 
sensitivity  matrices  were  developed  and  evaluated  to  minimize  the  matrix  condition 
number.  The  best-conditioned  sensitivity  matrix  was  developed  from  a  combination 
of  two-two  spring  sensitivity  matrices.  Similarly,  each  simulation  developed  a  well- 
conditioned  composite  sensitivity  matrix. 

Each  set  of  conditions  was  synthesized  into  FRF  data  gathered  from  the 
prototype  beam  the  finite  element  model  was  based  on.  The  resulting  element  wise 
correction  was  compared  to  a  previous  thesis  where  the  researcher  hand-picked 
modes  from  synthesized  pin  connections  based  on  the  QR  decomposition  of  the 
resulting  composite  sensitivity  matrix.  The  resulting  frequencies  are  displayed  in 
Table  22  through  Table  25  and  are  compared  to  previous  work  developed  in  [11], 
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Table  22.  Single  Mass  Solution  Comparison  to  Previous  Work 


Measured 
Freq  (Hz) 

Single  Spring 
Updated 
FEM  Freq 
(Hz) 

Single 

Spring 

Absolute 

Error 

Single 
Spring  % 
Relative 

Error 

Damanakis' 
Updated 
Frequencies 
(Hz)  [11] 

Pin 

Synthesis 

Absolute 

Error 

Pin 

Synthesis 
%  Relative 

Error 

Pin  vs 
Spring  % 
Relative 

Error 

13.8128 

13.996 

-0.1832 

1.3262 

13.8090 

-0.0038 

0.0275 

1.29870 

38.2736 

36.7333 

1.5403 

4.0244 

38.0604 

-0.2132 

0.5569 

3.4675 

75.0299 

74.5518 

0.4781 

0.6372 

74.6133 

-0.4166 

0.5553 

0.0819 

124.0133 

124.5348 

-0.5215 

0.4205 

123.3636 

-0.6497 

0.5239 

-0.1033 

185.2154 

186.5883 

-1.3729 

0.7412 

184.3768 

-0.8386 

0.4528 

0.2885 

258.5825 

261.7421 

-3.1596 

1.2219 

257.7584 

-0.8241 

0.3187 

0.9032 

344.2565 

348.7656 

-4.5091 

1.3098 

343.6795 

-0.5770 

0.1676 

1.1422 

442.2589 

448.5534 

-6.2945 

1.4233 

442.3850 

0.1261 

0.0285 

1.3948 

552.6715 

563.0675 

-10.3960 

1.8811 

554.1845 

1.5130 

0.2738 

1.6073 

674.5217 

687.8183 

-13.2966 

1.9713 

679.3893 

4.8676 

0.7216 

1.2496 

808.3990 

828.9641 

-20.5652 

2.5439 

818.0711 

9.6721 

1.1965 

1.3475 

953.9388 

980.0004 

-26.0616 

2.7320 

969.0386 

15.0998 

1.5829 

1.1491 

Table  23.  Single  Mass  Solution  Comparison  to  Previous  Work 


Measured 
Freq  (Hz) 

Single  Mass 
Updated 
FEM  Freq 
(Hz) 

Single 

Mass 

Absolute 

Error 

Single 
Mass  % 

Relative 

Error 

Damanakis' 
Updated 
Frequencies 
(Hz)  [11] 

Pin 

Synthesis 

Absolute 

Error 

Pin 

Synthesis 
%  Relative 

Error 

Pin  vs 

Mass  % 

Relative 

Error 

13.8128 

13.6115 

0.2013 

1.4574 

13.8090 

-0.0038 

0.0275 

1.42989 

38.2736 

35.7128 

2.5608 

6.6907 

38.0604 

-0.2132 

0.5569 

6.1338 

75.0299 

71.4144 

3.6155 

4.8188 

74.6133 

-0.4166 

0.5553 

4.2635 

124.0133 

118.0948 

5.9185 

4.7724 

123.3636 

-0.6497 

0.5239 

4.2486 

185.2154 

184.939 

0.2764 

0.1492 

184.3768 

-0.8386 

0.4528 

-0.3035 

258.5825 

257.2172 

1.3653 

0.5280 

257.7584 

-0.8241 

0.3187 

0.2093 

344.2565 

344.0446 

0.2119 

0.0616 

343.6795 

-0.5770 

0.1676 

-0.1061 

442.2589 

443.6277 

-1.3688 

0.3095 

442.3850 

0.1261 

0.0285 

0.2810 

552.6715 

553.5058 

-0.8343 

0.1510 

554.1845 

1.5130 

0.2738 

-0.1228 

674.5217 

679.6416 

-5.1199 

0.7590 

679.3893 

4.8676 

0.7216 

0.0374 

808.3990 

815.3904 

-6.9914 

0.8649 

818.0711 

9.6721 

1.1965 

-0.3316 

953.9388 

954.7796 

-0.8408 

0.0881 

969.0386 

15.0998 

1.5829 

-1.4947 
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Table  24.  Two  Spring  Solution  Comparison  to  Previous  Work 


Measured 
Freq  (Hz) 

Two 

Spring 

Updated 

FEM 

Freq  (Hz) 

Two 

Spring 

Absolute 

Error 

Two 

Spring 

% 

Relative 

Error 

Damanakis' 
Updated 
Frequencies 
(Hz)  [11] 

Pin 

Synthesis 

Absolute 

Error 

Pin 

Synthesis 
%  Relative 

Error 

Pin  vs 
Spring  % 
Relative 

Error 

13.8128 

13.7422 

0.0706 

0.5112 

13.8090 

-0.0038 

0.0275 

0.48366 

38.2736 

36.4948 

1.7788 

4.6475 

38.0604 

-0.2132 

0.5569 

4.0906 

75.0299 

74.8184 

0.2115 

0.2819 

74.6133 

-0.4166 

0.5553 

-0.2734 

124.0133 

120.4998 

3.5135 

2.8331 

123.3636 

-0.6497 

0.5239 

2.3093 

185.2154 

179.6513 

5.5641 

3.0041 

184.3768 

-0.8386 

0.4528 

2.5514 

258.5825 

249.6324 

8.9501 

3.4612 

257.7584 

-0.8241 

0.3187 

3.1425 

344.2565 

340.3857 

3.8708 

1.1244 

343.6795 

-0.5770 

0.1676 

0.9568 

442.2589 

443.0534 

-0.7945 

0.1797 

442.3850 

0.1261 

0.0285 

0.1511 

552.6715 

554.6007 

-1.9292 

0.3491 

554.1845 

1.5130 

0.2738 

0.0753 

674.5217 

679.7198 

-5.1981 

0.7706 

679.3893 

4.8676 

0.7216 

0.0490 

808.3990 

815.7007 

-7.3017 

0.9032 

818.0711 

9.6721 

1.1965 

-0.2932 

953.9388 

960.9866 

-7.0478 

0.7388 

969.0386 

15.0998 

1.5829 

-0.8441 

Table  25.  Two  Spring,  Low  Modes  Comparison  to  Previous  Work 


Measured 
Freq  (Hz) 

Two  Spring 
Low 

Updated 
FEM  Freq 
(Hz) 

Two 

Spring 

Low 

Absolute 

Error 

Two 
Spring 
Low  % 

Relative 

Error 

Damanakis' 
Updated 
Frequencies 
(Hz)  [11] 

Pin 

Synthesis 

Absolute 

Error 

Pin 

Synthesis 
%  Relative 

Error 

Pin  vs 
Spring  % 
Relative 

Error 

13.8128 

13.9022 

-0.0894 

0.6472 

13.8090 

-0.0038 

0.0275 

0.61962 

38.2736 

38.3917 

-0.1181 

0.3086 

38.0604 

-0.2132 

0.5569 

-0.2483 

75.0299 

75.1602 

-0.1303 

0.1736 

74.6133 

-0.4166 

0.5553 

-0.3817 

124.0133 

125.4511 

-1.4378 

1.1594 

123.3636 

-0.6497 

0.5239 

0.6356 

185.2154 

188.6115 

-3.3961 

1.8336 

184.3768 

-0.8386 

0.4528 

1.3808 

258.5825 

263.337 

-4.7545 

1.8387 

257.7584 

-0.8241 

0.3187 

1.5200 

344.2565 

351.5785 

-7.3220 

2.1269 

343.6795 

-0.5770 

0.1676 

1.9593 

442.2589 

454.3128 

-12.0539 

2.7255 

442.3850 

0.1261 

0.0285 

2.6970 

552.6715 

566.7458 

-14.0743 

2.5466 

554.1845 

1.5130 

0.2738 

2.2728 

674.5217 

694.1887 

-19.6670 

2.9157 

679.3893 

4.8676 

0.7216 

2.1941 

808.3990 

836.4128 

-28.0138 

3.4653 

818.0711 

9.6721 

1.1965 

2.2689 

953.9388 

991.8157 

-37.8769 

3.9706 

969.0386 

15.0998 

1.5829 

2.3877 
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As  can  be  seen  from  Table  22  through  Table  25,  the  synthesis  of  an 
artificial  pin  connection  and  time  intensive  mode  evaluation  still  provides  the  best 
model  updating;  however,  there  is  still  a  lot  of  information  that  can  be  gathered 
from  continuing  research  into  the  synthesis  of  an  artificial  spring  or  mass. 

B.  RECOMMENDATIONS 

Throughout  this  research,  any  spring  simulation  had  a  self-imposed  limit  of 
178.579  kg/mm  (10,000  Ib/in)  in  order  to  develop  an  understanding  of  the 
reactions  of  the  beam.  As  indicated  in  Figure  162  and  Figure  163,  there  may  be 
additional  information  for  the  sensitivity  matrix  with  spring  values  as  high  as 
2678.695  kg/mm  (150,000  Ib/in)  for  a  single  spring  system,  and  4464.491  kg/mm 
(250,000  Ib/in)  for  a  two  spring  system. 


Figure  1 62.  Effect  on  Sensitivity  Value  as  Single  Spring  Strength  Increases 
from  0-8928.984  kg/mm  (500,000  Ib/in)  on  Node  13 
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Figure  1 63.  Effect  on  Sensitivity  Value  as  Two  Spring  Strength  increases 
from  0-8928.984  kg/mm  (500,000  Ib/in)  on  Nodes  5  and  12 

equally 

As  indicated  in  Chapter  IV,  the  two-spring  synthesis  was  conducted  as  a 
proof  of  concept,  and  did  not  include  all  possibilities  for  evaluation.  There  are  105 
different  node  location  combinations,  and  multiple  spring  strengths  that  could 
produce  useable  FRF  and  sensitivity  matrices  from  each  set  of  nodes.  As  such, 
there  may  be  over  100,000  different  sensitivity  combinations  using  the  procedure 
detailed  here,  and  many  more  depending  on  how  many  modes  and  locations  are 
used  in  the  development  of  the  composite  sensitivity  matrix. 

Another  aspect  of  exploration  with  this  procedure  is  using  springs  of 
different  strengths  in  the  two  spring  synthesis.  As  indicated  in  Chapter  IV,  all  two- 
spring  synthesis  was  done  with  both  springs  of  equal  spring  strength.  By  varying 
spring  strength  separately,  the  dynamic  response  of  the  beam  should  alter 
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dramatically  and  provide  new  FRF  and  sensitivity  data  for  damage  detection  and 
model  updating. 

Ultimately,  because  there  is  so  much  information  available  from  this 
procedure,  a  program  should  be  developed  that  can  build  and  evaluate  the  rank 
and  condition  of  a  composite  sensitivity  matrix  based  on  either  existing  sensitivity 
matrices,  or  have  the  capability  to  build  new  matrices  as  required.  Such  a 
program  would  decrease  a  lot  of  the  time  required  in  developing  these  matrices 
and  reduce  the  subjectivity  in  development. 
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APPENDIX  A.  BEAM  MODIFICATION  BASE  CODE 


%%  Konze  -  Thesis 
clear, clc,clf, close  all 
format  shortG 

%  One  modification,  multiple  values,  calculates  S  for  each  condition  and 
%  pairs  down  based  on  selected  modes 

%%  Beam  Parameters 

%  Beam  Properties 
Beam.elem=16; 

Beam. ro=0. 0975/386. 4;  %lbm/inA3  (6061  Aluminum) 

Beam.E=10*10A6*.97;  %psi  with  0.9  correction  factor 

Beam.length=72;  %ft  to  in 

Beam.width=1 .93;  %in 

Beam. depth^O. 373;  %in 

Beam.dof=2;  %Degrees  of  freedom  per  node 

Beam.mass=Beam.length*Beam.depth*Beam.width*Beam.ro;  %lbm 

%Beam  as  described  by  Manos 

Beam.nodes=[0  0.5  5.57  10.64  15.71  20.78  25.85  30.92  35.99  41.06... 
46.1 3  51 .2  56.27  61 .34  66.41  71 .48  71 .98]; 


%  Beam. nodes=linspace(lelem_outer, Beam. length-lelem_outer, Beam. elem+1); 
BC=2;  %  Boundry  Condition.  1  =Cantilevered,  2=FreeFree 

Beam.l=Beam.width*Beam.depthA3/12;  %  Crossectional  Moment  of  Inertia 
Beam.gam=Beam.ro*Beam.depth*Beam.width;  %  Mass  per  unit  length 

%  Make  an  array  of  values  for  the  El  of  each  element.  Used  for  peturbance 
%  and  stiffness  matrix  building 
Eltmp=ones(1  ,Beam.elem); 

EI=(Beam.E*Beam.l); 

Beam.EI=EI*Eltmp; 


%  System  Limitations 

LFreq=1000;  %  Limiting  Frequency  (Hz) 


%  Artifical  Boundry  Condition  Characteristics 

ABC.SpringMag=500000;  %(lb/in)  Artifical  Spring  maximum  magnitude 
ABC. MassMag=0;%. 35/386. 4;  %(lbm)  Artifical  Mass  maximum  magnitude 

ABC.Node=[13];  %  Node  where  Artifical  Condition  is  placed 
ABC.dof=2.*ABC. Node-1 ;  %Specifies  the  translational  DOF  that  the  ABC  is  placed 

ABC.Interval=500;  %How  many  steps  to  go  from  1  to  the  magnitude  required 

if  ABC.SpringMag>0 

ABC. Spring  =  linspace(1  ,ABC.SpringMag,ABC.Interval);%  (Ib/ft) 

169 


else 

ABC.Spring=[]; 

end 

%  ABC.Spring=[1 0000]; 

%  ABC. lnterval=length(ABC. Spring); 

if  ABC.MassMag>0 

%  ABC. Mass  =  linspace(ABC. MassMag, 10/386.4, ABC. Interval);%  (Ib/ft) 
%  ABC. Mass  =  [.0001  .001  .01  .1  1 0  1 00  1 000  1 0000]/386.4  ; 
ABC.Mass=4/386.4; 
else 

ABC.Mass=[]; 

end 

ABC.Bmodes=0;  %  Baseline  higher  freq  modes  (+1 )  to  keep 
ABC.Mmodes=0;  %  Modified  higher  freq  modes  (+1 )  to  keep 


%  Damage 
Damage.  elem=5; 
Damage.  mag=5; 


%  Damaged  Element 
%  Damage  Magnitude  (percent) 


%%  Instrumentation  Properties 
%  Accelerometer 

total_accels=15;  %number  of  accelerometers  installed 
accel_weight=6. 5*0. 002204  ;%  gram  to  lbs 
h_accel=0.487  ;  %height  (in) 
d_accel=0.62  ;  %depth  (in) 

%Screw 

screw_weight=. 5*0. 002204  ;%  gram  to  lbs 
%Accelerometers 

lnst.accel_mass=accel_weight/386.4  ;  %lbs*sA2/in 

lnst.massl=lnst.accel_mass  /1 2  *(3*(d_accel/2)A2  +h_accelA2)  ;  %lbs  *  sA2  *in  Inertial  Mass 
%Screw 

screw_mass=screw_weight/386.4  ;  %lbs*sA2/in 
%Assy 

lnst.assy_mass=screw_mass+lnst.accel_mass  ;  %total  mass  of  screw-accel  assy 

lnst.nodes=[0.5  5.57  1 0.64  1 5.71  20.78  25.85  30.92  35.99  41 .06... 

46.1 3  51 .2  56.27  61 .34  66.41  71 .48]; 
for  r=1  :length(lnst. nodes) 

lnst.Node(r)=find(Beam.nodes==lnst.nodes(r)); 

end 

if  length(lnst.Node)~=length(lnst. nodes) 
display(‘need  to  fix  instrument/node  positions’) 
return 
end 

%  Remove  instrumentation  for  comparison  with  analytical  frequency  assement 
%  only 
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%  lnst.assy_mass=0; 
%  lnst.massl=0; 


%%  Frequencies  Measured 

Measured. freqs=[1 3.81 281 ,38.27358,75.02992,1 24.01 327,1 85.21 543,258.58254,344.25651 ,442. 
25887,552.67146,674.52174,808.39895,953.93877]; 


%%  Calculate  global  stiffness  and  mass  matricies 

[Beam.  KGIobal,  Beam.  MGIobal]=KandMbeam_Global(Beam,  Beam.  El,  Inst,  BC); 
%%  Analytical  frequencies  Assessment 

%  The  calculated  resonance  frequencies  are  compared  to  the  analytical 
%  solution  using  the  roots  to  the  chacteristic  equation  for  a  uniform, 

%  cantilevered  beam  as  described  using  the  formula  from  the  Fudamentls  of 
%  Structural  Dynamics,  2ed.,  example  13.3 

%  Analytic  Frequencies 
lamL(1)=1 .8751 ; 
laml_(2)=4.6941 ; 
laml_(3)=7.8548; 
laml_(4)=1 0.996; 


for  r=1  :length(laml_) 

om(r)=((lamL(r))A2/Beam.lengthA2)*(max(Beam.EI)/(Beam.ro*Beam.width*Beam.depth))A(1/2); 

end 

f_act=om/(2*pi()); 
if  BC==1 

KGIobal_BC=Beam. KGIobal; 

MGIobal_BC=Beam.MGIobal; 

else 

KGIobal_BC=Beam.KGIobal(Beam.dof+1  :end,Beam.dof+1  :end); 
MGIobal_BC=Beam.MGIobal(Beam.dof+1  :end,Beam.dof+1  :end); 
end 

[phi_BC,lam_BC,fn_BC]=eig_fcn(KGIobal_BC,MGIobal_BC); 

Analytical_err=(f_act’-fn_BC(1 :4))/f_act’*1 00; 

display(  Based  off  of  cantilevered  beam  analytical  assumption:') 

vars={‘Analytical_Freq’,’Calculated_Freq’,’Percent_Rel_Error’}; 

Comparison=table(f_act’,fn_BC(1 :4),Analytical_err(:, end), ’Variable’,  vars); 
display(Comparison) 


%%  Eigen  Values 
format  shortG 

[Beam.  phiFE,  Beam,  lam  FE,  Beam.  fnFE]=eig_fcn(Beam.  KGIobal,  Beam.  MGIobal); 


if  BC==1 

display(‘For  a  Cantilevered  beam,  the  calculated  resonance  frequencies  are:’) 
Freq_Calculated=Beam.fnFE’; 
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display(Beam.fnFE) 
elseif  BC==2 

%comaprison  against  measured  frequencies  is  only  vaid  for  a  free  free 
%beam 

display(‘For  a  free-free  beam,  with  ridged  body  modes  removed,  the  calculated  resonance 
frequencies  are:’) 

Freq_Calculated=Beam.fnFE(3:14)’; 

display(Beam.fnFE(3:end)) 

meas_err=abs((Measured.freqs’-Freq_Calculated’)/Measured.freqs’)*100; 

vars1={‘Measured_Freq’,’Calculated_Freq’,’Percent_Rel_Error’}; 

table(Measured.freqs’,Freq_Calculated’,meas_err(:,end),’Variable’,vars1) 

else 

display(‘Fix  the  BCs’) 
end 


format  short 


%%  System  Frequency  Limitations 

[Beam.  LimitedPhi,  Beam. Limited  Lam, Beam.  LimitedFn]=sys_limit(Beam.phiFE, Beam,  lam  FE,Bea 
m.fnFE,LFreq); 

%%  Sensitivity  Calculation 

[Beam.  S]=Sensitivity_fcn(Beam,  Inst,  BC); 


[Check]=Sensitivity_Check(Beam,lnst,BC); 
display(‘Sensitivity  Relative  Error  (in  percent)') 

Check 

display(‘The  largest  sensitivity  error  is  (percent)’) 
max(max(Check)) 

%Beam.S=Beam.S(3:1 6,:); 

R_C.Base.rank=rank(Beam.S); 

R_C.Base.cond=cond(Beam.S); 

%%  Apply  simulated  spring  and  mass  for  FRF  analysis 
[Delta_fn]=ABC_FRF(Beam,ABC,BC,LFreq); 

%%  Apply  simulated  spring  and  mass  for  Damage  Detection 
Pivot.  CondS=ones(1, 14); 

% 

%  cnt=1 
cnt=0; 

for  CntElem=2:15 
cnt=cnt+1 ; 

Damage.  elem=CntElem; 

fprintf(  For  element  %g  \n\n’, Damage. elem) 

if  LFreq>=Beam.fnFE(end) 
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[phiTest,lamTest,fnTest]=damage_beam(Beam,lnst,BC,  Damage); 
[L_phitest,L_lamtest,L_fntest]=sysJimit(phiTest,lamTest,fnTest,  LFreq); 
dlam=L_lamtest-Beam.LimitedLam; 


o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o 

%dlam=dlam(3:16); 

o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o 


%  deltaP=inv(Beam.S)*dlam; 

deltaP=Beam.S\dlam; 

%  deltaP=Beam.S\dlam(1 :14); 

figure 

TrueDeltaP=Beam.EI’; 

TrueDeltaP(Damage.elem)=TrueDeltaP(Damage.elem)+TrueDeltaP(Damage.elem)*Damage.ma 

g/100; 

T  rueDeltaP=Beam.EI’-T  rueDeltaP; 

TrueDelta_P=TrueDeltaP(2:end-1); 

%  Damage  Detection  Plots 
barThis  =  [TrueDelta_P  deltaP]; 
bar(barThis) 

title(' Plot:  Model  Errors  vs  True  Errors  Using  All  Modes’) 
legend(‘True’, ’Indicated’, ’Location’, ’Best) 
xlabel(‘Element  #’) 

%axis([0  18-10000  10000]) 
else 

%  Perform  Damage  Detection  with  Spring 
tic 

[S_s,dlam_s,  Pivot]=ABC_Sensistivity(Beam,  ABC,  LFreq,  Inst,  Damage,  BC,  Pivot); 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

%For  testing  a  composite  spring  system  and  using  excel 
%  [S_s,dlam_s,Pivot]=ABC_CompositeS(Beam,  ABC,  LFreq,  Inst,  Damage,  BC,  Pivot); 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


deltaP_s=S_s\dlam_s; 

R_C.spring_rank(cnt)=rank(S_s); 

R_C.spring_cond(cnt)=cond(S_s); 

o/  o/  o/ 

/o  /o  /o 

SpringRank=rank(S_s) 

Pivot.CondS(cnt)=Pivot.CondS(cnt)*cond(S_s); 

SpringTime=toc 

o/  o/  o/ 

/o  /o  /o 

%Perform  Damage  Detection  with  Pin 
%  tic 

%  [S_p,dlam_p,  Pivot,  pin,  Pivot,  pin  J]=ABC_Pinned(Beam,  LFreq,  Inst,  Damage,  BC,  ABC); 

%  deltaP_p=S_p\dlam_p; 

%  R_C.pin_rank(cnt)=rank(S_p); 

%  R_C.pin_cond(cnt)=cond(S_p); 
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% 


PinTime=toc 


figure 

TrueDeltaP=Beam.EI’; 

TrueDeltaP(Damage.elem)=TrueDeltaP(Damage.elem)+TrueDeltaP(Damage.elem)*Damage.ma 

g/100; 

TrueDeltaP=Beam.Er-TrueDeltaP; 

T  rueDelta_P=T  rueDeltaP(2:end-1 ); 

%  Damage  Detection  Plots 
%  barThis  =  [TrueDelta_P  deltaP_s  deltaP_p]; 
barThis  =  [TrueDelta_P  deltaP_s]; 
bar(barThis) 

title(' Plot:  Model  Errors  vs  True  Errors  with  Limited  Freq  of  1000Hz’) 
legend(‘True’, ’Spring’, ’Pin’, ’Location’, ’Best’) 
xlabel(‘Element  #’) 

%axis([0  18-10000  10000]) 
end 
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APPENDIX  B.  STIFFNESS  AND  MASS  MATRIX  BUILDING 

FUNCTION 


function  [KGIobal,MGIobal]=KandMbeam_Global(Beam,EI,lnst,BC) 


% 

%  Purpose  :  To  develop  a  Global  stiffness  and  mass  matrix  for  a  beam 
%  Method  (Two  parts)  :  (Part  1)  Using  the  shape  functions  for  a  beam, 

%  perform  symbolic  integration  and  differentiation  to  develop  the  stiffness 
%  and  mass  matrix  for  the  element(Longer  processing  time). 

%  (Part  2)  Using  known  solutions  to  the  integrations,  develop  the  stiffness 
%  and  mass  matrix  for  the  element  (faster  processing  time) 

% 

%  Variable  Description: 

%  INPUT  parameters: 

%  Beam:  Properties  of  the  beam  with  the  following  entries: 

%  length 

%  width 

%  depth 

%  elem:  Number  of  elements  in  the  beam 

%  nodes:  Array  of  nodal  locations 

%  El:  Modulus  of  Elasticity  multiplied  by  the  crossectional 
%  moment  of  inertia 

%  Inst:  Properties  of  the  instruments  with  the  following  entries: 

%  assy_mass:  Any  additional  mass  located  at  the  nodes  of  the  beam 
%  (e.g.  Acceleromerters) 

%  massl:  Innertal  mass  of  the  additional  mass,  againlocated  at 
%  the  nodes 

% 

%  OUTPUT  parameters: 

%  KGIobal:  Global  Stiffness  Matrix,  no  boundry  conditions 

%  MGIobal:  Global  Mass  Matrix,  no  boundry  conditions 

% - 

syms  x  y  z 


%  %%  Part  1  (uncomment  to  use) 

%  sf(1  )=1  -3*(x/lelem)A2+2*(x/lelem)A3; 

%  sf(2)=x*(1  -2*(x/lelem)+(x/lelem)A2); 

%  sf(3)=3*(x/lelem)A2-2*(x/lelem)A3; 

%  sf(4)=-lelem*(x/lelem)A2+lelem*(x/lelem)A3; 

% 

% 

%  for  i=1  :length(sf) 

%  for  j=1  :length(sf) 

%  meiem(i,j)=gam*int(sf(i)*sf(j),x,0,ielem); 

%  keiem(i,j)=EI(j)*int(diff(diff(sf(i),x),x)*diff(diff(sf(j),x),x),x,0,ielem); 

%  end 
%  end 
% 

%  r=1 ;  %element  starting  counter 
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%  v=4;  %element  ending  counter 

% 

%  KGIobal=zeros((elem*dof)+dof);  %KGIobal  preallocation 
%  MGIobal=zeros((elem*dof)+dof);  %MGIobal  preallocation 

% 

%  while  r<=elem*dof 

% 

% 

%  %  Assemble  Stiffness  and  Mass  Global  matricies 

%  KGIobal(r:v,r:v)=KGIobal(r:v,r:v)+kelem; 

%  MGIobal(r:v,r:v)=MGIobal(r:v,r:v)+melem; 

% 

%  %  Moving  the  number  of  degrees  of  freedom  for  per  node  in  rows  and 

%  %  columns  after  each  loop 

%  r=r+dof; 

%  v=v+dof; 

%  n=n+1 ; 

%  end 


%%%  Part  2  (uncomment  to  use) 

%establish  counters  and  preallocate  global  matricies 
n=1 ;  %  Node  starting  counter 

r=1 ;  %element  starting  counter 
v=4;  %element  ending  counter 

KGIobal=zeros((Beam.elem*Beam.dof)+Beam.dof);  %KGIobal  preallocation 
MGIobal=zeros((Beam.elem*Beam.dof)+Beam.dof);  %MGIobal  preallocation 

while  v<=length(Beam.nodes)*Beam.dof 

lelem=Beam.nodes(n+1)-Beam.nodes(n); 

kelem=(EI(n)/(lelemA3))*... 

[12  6*lelem  -12  6*lelem; 

6*lelem  4*lelemA2  -6*lelem  2*lelemA2; 

-12  -6*lelem  12  -6*lelem; 

6*lelem  2*lelemA2  -6*lelem  4*lelemA2]; 

melem=(Beam.ro*lelem*Beam.width*Beam. depth/420)*... 

[  156  22*lelem  54  -13*lelem; 

22*lelem  4*lelemA2  13*lelem  -3*lelemA2; 

54  13*lelem  156  -22*lelem; 

-13*lelem  -3*lelemA2  -22*lelem  4*lelemA2[; 


%  Assemble  Stiffness  and  Mass  Global  matricies 

KGIobal(r:v,r:v)=KGIobal(r:v,r:v)+kelem; 

MGIobal(r:v,r:v)=MGIobal(r:v,r:v)+melem; 

%  Moving  the  number  of  degrees  of  freedom  for  per  node  in  rows  and 

%  columns  after  each  loop 

r=r+Beam.dof; 

v=v+Beam.dof; 

n=n+1 ; 
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end 


%  Add  point  and  inertial  loads  due  to  accelerometers  (no  accelerometer  at 
%  the  ends  of  the  beam) 
for  rr=1  :length(lnst. nodes); 

l_node=find(lnst.nodes(rr)==Beam.  nodes); 
dof1=l_node*2-1 ; 
dof2=dof1+1 ; 

MGIobal(dof1  ,dof1  )=MGIobal(dof1  ,dof1  )+lnst.assy_mass; 
MGIobal(dof2,dof2)=MGIobal(dof2,dof2)+lnst.massl; 

end 


if  BC==1 

KGIobal=KGIobal(3:end,3:end); 

MGIobal=MGIobal(3:end,3:end); 
%  nodes=nodes(find(nodes>0)); 
elseif  BC==2 

KGIobal=KGIobal; 

MGIobal=MGIobal; 

%  nodes=nodes; 
else 

display(‘Fix  the  BCs’) 
end 
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APPENDIX  C.  EIGENVALUE  FUNCTION 


function  [phi,lam,fn]=eig_fcn(KGIobal,MGIobal) 


o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


%  Purpose: 

%  Take  the  input  stiffness  and  mass  matricies  and  output  the  eignenvectors, 
%  eigenvalues,  and  frequencies  of  resonance. 

%  INPUT  parameters: 

%  KGIobal:  Global  Stiffness  Matrix 

%  MGIobal:  Global  Mass  Matrix 

%  OUTPUT  parameters: 

%  phi:  Eigenvectors  of  the  matricies 

%  lam:  Eigenvalues  of  the  matricies 

%  fn:  Resonance  frequencies  of  the  system  (in  Hz) 


o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


[phi, lam]  =  eig(KGIobal, MGIobal); 
lam  =  diag(lam); 
wn  =  real(sqrt(lam)); 
fn  =  wn/(2*pi()); 


179 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


180 


APPENDIX  D.  SYSTEM  FREQUENCY  LIMITING  FUNCTION 


function  [L_phi,LJam,L_fn]=sysJimit(phi,lam,fn,LFreq) 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


%  Purpose: 

%  output  an  updated  eigenvalue  and  eigenvector  based  on  the  limiting 
%  frequency.  The  limiting  frequency  may  be  due  to  system  bandwidth,  and 
%  will  affect  the  readable  eigenvalues  and  eigenvectors.  It  can  be  assumed 
%  that  modes  and  resonance  frequencies  above  this  limiting  frequency  are 
%  nolonger  “viewable,”  and  may  be  unuseable  in  future  calculations. 

% 

%  This  function  also  removes  free  body  modes,  should  they  exist. 

% 


%  Input  parameters: 

%  phi:  Eigenvectors  of  the  matricies 

%  lam:  Eigenvalues  of  the  matricies 

%  fn:  Resonance  frequencied  of  the  system  (in  Hz) 
%  LFreq:  Limiting  frequency 

% 


%  Output  Parameters 
%  L  phi:  “viewable”  Eigenvectors 
%  LJam:  “viewable”  Eigenvalues  of  the  matricies 
%  L_fn:  “viewable"  Resonance  frequencies  of  the  system  (in  Hz) 


o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


%Below  loop  limits  the  usabel  modes  and  frequencies  based  off  of  LFreq 
r=1 ; 

if  LFreq>=fn(end) 

L_phi=phi; 

L_iam=lam; 

L_fn=fn; 

else 

while  fn(r)<=LFreq  &&  r<=length(fn) 

L_phi(:,r)=phi(:,r); 

L_lam(r,1)=lam(r); 

L_fn(r,1)=fn(r); 
r=r+1 ; 
end 
end 

%  Below  loop  removes  free  body  modes, 
for  rr=1 :4 

if  L_fn(1  )<=1 0A(-2) 

L_phi=L_phi(:,2:end); 

L_lam=LJam(2:end); 

L_fn=L_fn(2:end); 

end 

end 
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APPENDIX  E.  SENSITIVITY  MATRIX  BUILDING  FUNCTION 


function  [S]=Sensitivity_fcn(Beam, Inst, BC, ABC, cnt,LFreq, node) 

% 

%  Purpose  :  To  develop  an  Eigen  Sensitivity  Matrix  based  on  changes  to  the 
%  El  values.  The  Sensitivity  matrix  assists  in  predicting  changes  to 
%  resonance  frequencies  based  on  changes  to  specific  parameters  in  the 
%  system.  This  maxrix  is  developed  by  temporarily  perturbing  the  system  and 
%  recording  how  the  finite  element  model  reacts. 

% 

%  Variable  Description: 

%  INPUT  parameters: 

%  Beam:  Properties  of  the  beam  with  the  following  entries: 

%  phiFE  or  LimitedPhi:  Eigenvectors  of  the  Stiffness  and  Mass  matricies 
%  KGIobal:  Global  Stiffness  Matrix 
%  El:  Modulus  of  Elasticity  multiplied  by  the  crossectional 
%  moment  of  inertia 

%  length  :  length  of  the  beam 
%  width  :  Width  of  the  beam 
%  depth  :  Depth  of  the  beam 
%  elem:  Number  of  elements  in  the  beam 
%  nodes:  Array  of  nodal  locations 

%  Inst:  Properties  of  the  instruments  with  te  following  entries: 

%  assy_mass:  Any  additional  mass  located  at  the  nodes  of  the  beam 
%  (e.g.  Acceleromerters) 

%  massl:  Innertal  mass  of  the  additional  mass,  againlocated  at 
%  the  nodes 

%  Node:  (optional)  Node  that  a  simulated  pin  connection  occurs  on. 

%  OUTPUT  parameters: 

%  S:  Eigen  Sensitivity  matrix 

% 

p=exist('node’); 

c=exist(ABC’); 

PercentPerturb  =  1 ;  %  finite  difference  step  size  for  calculating  sensitivity  (unrelated  to  actual 
error  in  FEM). 

NumModes=min(size(Beam.  LimitedPhi)); 

%  NumModes=14; 

NumSprings=Beam.elem; 

%  S  =  zeros(NumModes,NumSprings);  %  initialize  S  matrix 

for  iCntElem  =  1  :  NumSprings 

%  KGIobal_Tmp  =  Beam. KGIobal; 

%Develop  an  array  for  the  peturbed  bar  using  the  El 
EI_pet=Beam.EI; 

EI_tmp=EI_pet(iCntElem)*(PercentPerturb/100); 

EI_pet(iCntElem)=EI_pet(iCntElem)+EI_tmp; 


[KGIobal_Tmp,MGIobal_Tmp]=KandMbeam_Global(Beam,EI_pet,lnst,BC); 
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if  c==1 


if  ABC.MassMag>0 
mdelta=ABC.Mass(cnt); 

MGIobal_Tmp(ABC.dof,ABC.dof)=MGIobal_Tmp(ABC.dof,ABC.dof)+mdelta; 

end 


if  ABC.SpringMag>0 
kdelta=ABC.Spring(cnt); 

KGIobal_Tmp(ABC.dof,ABC.dof)=KGIobal_Tmp(ABC.dof,ABC.dof)+kdelta; 

end 

end 


if  p==1 


[KGIobal_Tmp,MGIobal_Tmp]=Pinned_Dof(KGIobal_Tmp,MGIobal_Tmp,node); 
[KGIobal_Pin,  MGIobal_Pin]=Pinned_Dof(Beam.  KGIobal,  Beam.  MGIobal,  node); 

[phi_mod,lam_mod,fn_mod]=eig_fcn(KGIobal_Pin,MGIobal_Pin); 

[L_phiPin,L_lamPin,L_fnPin]=sys_limit(phi_mod,lam_mod,fn_mod,LFreq); 

dKdk=(KGIobal_Tmp-KGIobal_Pin)/EI_tmp; 

NumModes=length(LJamPin); 

for  iCntModes  =  1  :  NumModes  %  Using  all  modes  of  the  system  (NumModes  =  nP  =  #rows 
of  [S]) 

S(iCntModes,iCntElem)  =  L_p h i P i n ( : , iCntModes)’  *  dKdk  *  L_p h i P i n ( : , iCntModes); 
end 
else 

dKdk=(KGIobal_Tmp-Beam.KGIobal)/EI_tmp; 

for  iCntModes  =  1  :  NumModes  %  Using  all  modes  of  the  system  (NumModes  =  nP  =  #rows 
of  [S]) 

S(iCntModes,iCntElem)  =  Beam. LimitedPhi(:, iCntModes)’  *  dKdk  * 

Beam.  LimitedPhi(:,  iCntModes); 
end 
end 

end 

S=S(:,2:end-1); 
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APPENDIX  F.  SENSITIVITY  MATRIX  CHECK  FUNCTION 


function  [RelErr]=Sensitivity_Check(Beam,lnst,BC) 

% 

%  Purpose  :  To  check  the  accuracy  of  an  Eigen  Sensitivity  Matrix  based  on 
%  changes  to  the  El  values.  The  Sensitivity  matrix  assists  in  predicting 
%  changes  to  resonance  frequencies  based  on  changes  to  specific  parameters 
%  in  the  system.  The  check  is  accomplished  by  comparing  the  difference 
%  between  the  the  resonance  frequency  prediction  to  change  using  the  Eigen 
%  Sensitivity  Matrix  to  a  system  with  the  change  installed. 

% 

%  Variable  Description: 

%  INPUT  parameters: 

%  Beam:  Properties  of  the  beam  with  the  following  entries: 

%  length 

%  width 

%  depth 

%  ro:  material  density 

%  elem:  Number  of  elements  in  the  beam 
%  S:  Eigen  Sensitivity  matrix 

%  nodes:  Array  of  nodal  locations 

%  El:  Modulus  of  Elasticity  multiplied  by  the  crossectional 
%  moment  of  inertia 

%  lam_base:Eigenvectors  of  the  Stiffness  and  Mass  matricies 
%  KGIobal:  Global  Stiffness  Matrix 

%  MGIobal:  Global  Mass  Matrix 

%  Inst:  Proterties  of  the  instrumentation  with  the  following  entries: 

%  assy_mass:  Any  additional  mass  located  at  the  nodes  of  the  beam 
%  (e.g.  Acceleromerters) 

%  massl:  Innertal  mass  of  the  additional  mass,  againlocated  at 
%  the  nodes 

%  BC:  Boundry  Conditions  enforced  (cantilevared  or  free  free) 

%  OUTPUT  parameters: 

%  RelErr:  Relative  Error  of  the  predicted  eigenvalue  compared  to 
%  the  eigen  value  with  the  temporary  modification  installed 


PercentChange  =  1 ;  %  finite  difference  step  size  for  calculating  sensitivity  (unrelated  to  actual 
error  in  FEM). 

NumModes=length(Beam.  Limited  Lam); 

%  NumModes=14; 

NumSprings=Beam.elem-2; 


RelErr  =  zeros(NumModes,NumSprings);  %  initialize  error  check  matrix 
for  iCntElem  =  1  :  NumSprings 

KGIobal_Tmp  =  Beam. KGIobal; 

EI_pet=Beam.EI; 

EI_tmp=Beam.EI(iCntElem)*(PercentChange/100); 
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EI_pet(iCntElem)=EI_pet(iCntElem)+EI_tmp; 
[KGIobal_Tmp,~]=KandMbeam_Global(Beam,EI_pet,lnst,BC); 
%  kcheck=KGIobal_Tmp-KGIobal; 


[~,LamChanged,~]=eig_fcn(KGIobal_Tmp,Beam.MGIobal); 

%Remove  Free  Body  Modes 
if  LamChanged(1)<=10A(-2) 

lamChanged=LamChanged(3:end); 

else 

lamChanged=LamChanged; 

end 

for  iCntModes  =  1  :  NumModes 

%  Predict  the  eigenvalue  (and  the  frequency  by  extension  )  of  the 
%  system  using  the  base  system  eigenvalues,  the  sensitivity  matrix, 

%  and  the  total  change. 

%  lamPredict=lam_base(iCntModes)+S(iCntModes,iCntElem)*EI_pet(iCntElem); 

lamPredict=Beam.LimitedLam(iCntModes)+Beam.S(iCntModes,iCntElem)*EI_tmp; 

%Ftelative  Error  of  the  sensitivity  matrix  in  percent 
RelErr(iCntModes,iCntElem)=abs((lamChanged(iCntModes)- 
lamPredict)/lamChanged(iCntModes))*100; 

end 

end 
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APPENDIX  G.  FRF  SYNTHESIS  FUNCTION 


function  [fn_mod_record,Beam]=ABC_FF!F(Beam,ABC,BC,LFreq) 


o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


%  Variable  Description: 

%  INPUT  parameters: 

%  Beam:  Properties  of  the  beam  with  the  following  entries: 

%  phiFE  or  LimitedPhi:  Eigenvectors  of  the  Stiffness  and  Mass  matricies 
%  KGIobal:  Global  Stiffness  Matrix 

%  El:  Modulus  of  Elasticity  multiplied  by  the  crossectional 
%  moment  of  inertia 

%  length  :  length  of  the  beam 
%  width  :  Width  of  the  beam 
%  depth  :  Depth  of  the  beam 
%  elem:  Number  of  elements  in  the  beam 
%  nodes:  Array  of  nodal  locations 

%  ABC:  Properties  of  the  Artifical  Boundary  Conditions 
%  SpringMag:  The  magnitude  of  artifical  spring-to-ground  stiffness 
%  MassMag:  The  magnitude  of  artifical  dimensionless  mass 
%  Node:  Node  to  apply  the  Artifical  beam  modifications 
%  dof:Degree  of  Freedom  associated  to  the  Node  for  Beam 
%  Modifications 

%  Interval:  How  many  steps  from  1  to  the  magnitude  of 
%  modification 

%  Spring:  An  array  of  spring  values  going  from  1  to  the  magnitude 
%  of  modification  with  ABC. Interval  steps 

%  Mass:  An  array  of  mass  values  going  from  1  to  the  magnitude 
%  of  modification  with  ABC. Interval  steps 

%  Inst:  Properties  of  the  instruments  with  te  following  entries: 

%  assy_mass:  Any  additional  mass  located  at  the  nodes  of  the  beam 
%  (e.g.  Acceleromerters) 

%  massl:  Innertal  mass  of  the  additional  mass,  againlocated  at 
%  the  nodes 

%  Damage:  Properties  of  the  damage  to  the  beam 
%  elem:  Element  that  damage  occurs  on 
%  Mag:  Manitude  of  damage 
%  BC:  Boundary  Condition  applied  to  the  beam 

%  LFreq:  Limiting  frequency  of  the  system.  Modes  and  resonances  above 
%  this  frequency  are  not  included  in  calculation 

% 


%  OUTPUT  parameters: 

%  fn_mod_record:  A  record  of  the  limited  frequencies  of  the  modified  beam 
%  for  every  step  of  ABC  interval, including  the  baseline 

%  beam 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


%%  Structural  Synthisis  Properties  for  Artifical  Boundry  Condition 


Freqs  =  1  :  0.1:  1000; 

fn_mod=zeros(length(Beam.LimitedFn),  ABC.  Interval); 
fn_record=zeros(length(Beam.LimitedFn)+1  ,ABC.Interval+1); 


if  length(ABC.Node)>1 
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ABC.Node=ABC.Node(1 ); 
ABC.dof=ABC.dof(1 ); 
end 


if  BC==1 

nodes=Beam.nodes(2:end); 

cset=ABC.dof; 

iset=[0:cset-1  ,cset+1 :2*length(nodes)]; 
iset=iset(2:end); 
aset=[cset  iset]; 
else 

nodes=Beam. nodes; 
cset=ABC.dof; 

iset=[0:cset-1  ,cset+1 :2*length(nodes)]; 
iset=iset(2:end); 
aset=[cset  iset]; 
end 

%%  Inital  FRF  Data 
cnt=0; 

for  freq  =  Freqs 

cnt  =  cnt  +  1 ;  %  counter  for  plot  vectors 
omega  =  2*pi*freq;  %  forcing  frequency  in  rad/sec 

Z  =  Beam.KGIobal-omegaA2*Beam.MGIobal;  %  impedance  matrix  at  frequency  omega 
Beam. cnt)  =  inv(Z);  %  FRF  matrix  at  frequency  omega 
H1 1  inital(cnt)  =  Beam.H(2,2,cnt); 

Hcs_inital(cnt)=Beam.H(cset,cset,cnt); 

His_inital(cnt)=Beam.H(iset(end-1),iset(end-1),cnt); 

end 

%%  Pin  ABC  FRF  Synthesis  data 


cnt=0; 

for  freq  =  Freqs 
cnt  =  cnt  +  1 ; 

Hii  =  Beam. H(iset, iset, cnt); 

Hie  =  Beam.  H(iset,cset,  cnt); 

Hci  =  Beam. H(cset, iset, cnt); 

Hcc  =  Beam. H(cset,cset, cnt); 

%  HiiStar  =  Hii  -  Hie  *  (Hcc\Hci); 

HStar_pin(:,:,cnt)  =  Hii  -  Hie  *  inv(Hcc)  *Hci; 

HiiStarpin(cnt)  =  HStar_pin(iset(end-1),iset(end-1),cnt); 

HccStar_pin(cnt)  =  HStar_pin(cset,cset,cnt); 
end 

%%  Pin  FRF  Check 

[KGIobal_Pin, MGIobal_Pin]=Pinned_Dof(Beam.  KGIobal,  Beam.  MGIobal,  ABC. Node); 


cnt=0; 

for  freq  =  Freqs 

cnt  =  cnt  +  1 ;  %  counter  for  plot  vectors 
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omega  =  2*pi*freq;  %  forcing  frequency  in  rad/sec 

Z  =  KGIobal_Pin-omegaA2*MGIobal_Pin;  %  impedance  matrix  at  frequency  omega 
H_pin(:,:,cnt)  =  inv(Z);  %  FRF  matrix  at  frequency  omega 
HI  1_pin_chk(cnt)  =  H_pin(1,1,cnt); 

Hcs_pin_chk(cnt)=H_pin(cset,cset,cnt); 

His_pin_chk(cnt)=H_pin(iset(end)-1  ,iset(end)-1  ,cnt); 
end 


%%  Plots 

figure(2) 
subplot(21 1) 

suptitle(  Inital  vs  Pin  FRF  response’) 

plot(Freqs,20*log10(abs(squeeze(Beam.H(iset(end)-1  ,iset(end)-1, :))))) 
hold  on 

plot(Freqs,20*log10(abs(HiiStar_pin))) 

axis([0  150-160  40]) 

legend(' Inital’, ’Pinned’, ’Location’, ’Best') 

title([‘|H_{‘,num2str(iset(end)-1 ),’  ‘,num2str(iset(end)-1),’}(\Omega)|’]) 
xlabel(‘Frequency  (Hz)’) 
subplot(21 2) 

plot(Freqs,angle(squeeze(Beam.H(iset(end)-1  ,iset(end)-1  ,:)))*1 80/pi) 
hold  on 

plot(Freqs,angle(HiiStar_pin)*1 80/pi) 
title('Phase  \Theta_{1 1 }  (\Omega)’) 

figure(3) 
subplot(21 1) 

suptitle(  Pin  FRF  response  check’) 
plot(Freqs,20*log10(abs(HccStar_pin))) 
hold  on 

plot(Freqs,20*log10(abs(Hcs_pin_chk)),’-’) 
legend('Check’,  ’Synthesized’,  ’Location’,  ’Best') 

title([‘ |H_{‘,num2str(iset(end)-1 ),’  ‘,num2str(iset(end)-1),’}(\Omega)|for  pinned  DOF  Check’]) 
xlabel(‘Frequency  (Hz)’) 
subplot(21 2) 

plot(Freqs,angle(HccStar_pin)*1 80/pi) 
hold  on 

plot(Freqs,angle(Hcs_pin_chk)*1 80/pi,’-') 
title('Phase  \Theta_{is  is]  (\Omega)’) 


%% 

fig_num=5; 


for  iCntlter=1  :ABC. Interval 

°/o°/o  ABC  Mass 

%  Add  point  loads  due  to  Artifical  Mass 
MGIobal_Mod=Beam.MGIobal; 
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if  ABC.Mass>0 
mdelta=ABC.Mass(iCntlter); 

MGIobal_Mod(ABC.dof,ABC.dof)=Beam.MGIobal(ABC.dof,ABC.dof)+mdelta; 

end 

%%  ABC  Spring 
KGIobal_Mod=Beam.KGIobal; 
if  ABC.Spring>0 

kdelta=ABC.Spring(iCntlter); 

KGIobal_Mod(ABC.dof,ABC.dof)=KGIobal_Mod(ABC.dof,ABC.dof)+kdelta; 

end 

%%  Beam  Modification  Integration  Sythisis  Check 
[phi,lam,fn]=eig_fcn(KGIobal_Mod,MGIobal_Mod); 
[L_phitest,L_lamtest,L_fntest]=sys_lirnit(phi,lam,fn,LFreq); 


fn_mod(1  :length(L_fntest),iCntlter)=L_fntest; 

fn_mod_record(1  :length(Beam.LimitedFn),iCntlter)=fn(1  :length(Beam.LimitedFn)); 
if  ABC.Spring~=0 

vars(iCntlter+1)={sprintf('%s%s%s’,  ’Spring’,  num2str(ABC.Spring(iCntlter)),’lb/in’)}; 
else 

vars(iCntlter+1)={sprintf('%s%s%s’,  'Mass’,  num2str(ABC.Mass(iCntlter)*386.4),  ’lb’)}; 
end 


cnt  =  0; 

Flmod=zeros(length(nodes)*Beam.dof,length(nodes)*Beam.dof,length(Freqs)); 

Hcscheck=zeros(1,length(Freqs)); 

Hischeck=zeros(1,length(Freqs)); 
for  freq  =  Freqs 

cnt  =  cnt  +  1 ;  %  counter  for  plot  vectors 

omega  =  2*pi*freq;  %  forcing  frequency  in  rad/sec 

Z  =  KGIobal_Mod-(omegaA2*MGIobal_Mod);  %  Inversion  matrix  at  frequency  omega 
Hmod(:,:,cnt)  =  inv(Z); 

Hcscheck(cnt)=Hmod(cset,cset,cnt); 

Hischeck(cnt)=Hmod(3,3,cnt); 

end 


%%  Beam  Modification  Integration  Synthisis 
cnt  =  0; 

Hstarcs=zeros(1  Jength(Freqs)); 
Hstaris=zeros(1  Jength(Freqs)); 
for  freq  =  Freqs 
cnt=cnt+1 ; 

Haa  =  Beam. H(aset,aset, cnt); 

Hac  =  Beam. H(aset,cset, cnt); 

Hca  =  Beam. H(cset,aset, cnt); 

Hcc  =  Beam. H(cset,cset, cnt); 

if  ABC.Spring~=0 
delZ=ABC.Spring(iCntlter); 
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else 

omega  =  2*pi*freq; 
delZ=-omegaA2*ABC.Mass(iCntlter); 
end 


Hstar(aset,aset,cnt)=Haa-Hac*inv(Hcc+(inv(delZ)))*Hca; 

Hstar_cssynth(cnt)=Hstar(cset,cset,cnt); 

Hstar_issynth(cnt)=Hstar(3,3,cnt); 


end 

%%  Plots 
figure(4) 

%  subplot(211) 

plot(Freqs,20*log10(abs(Hstar_cssynth))) 
hold  on 


% 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

%  %  Testing  only 

axis([0  150-160  40]) 
grid  on 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

%  title([‘ |H_{‘,num2str(iset(1 )),’  ‘,num2str(iset(1)),’}(\Omega)|’]) 

%  xlabel(‘Frequency  (Hz)’) 

%  subplot(212) 

%  plot(Freqs,angle(Hstar_cssynth)*1 80/pi) 

%  hold  on 

%  %  title(‘Phase  \Theta_{is  is}  (\Omega)’) 
figure 

subplot(211) 

suptitle(‘Modified  CSET  FRF  response  check’) 
plot(Freqs,20*log10(abs(Hcscheck)),’-') 
hold  on 

plot(Freqs,20*log10(abs(Hstar_cssynth)),’-.') 
legend(‘Check’, ’Synthesized’, ’Location’, ’Best’) 
if  ABC.Spring~=0 

title([‘|H_{‘,num2str(cset),’  ‘,num2str(cset),... 

‘}(\Omega)|for  k=‘,num2str(ABC.Spring(iCntlter)), ’in/lb  spring’]) 

else 

title([‘|H_{‘,num2str(cset),’  ‘,num2str(cset),... 

‘}(\Omega)|for  m=‘,num2str(ABC.Mass(iCntlter)),’lbm  mass’]) 
end 

xlabel('Frequency  (Hz)’) 
subplot(212) 

plot(Freqs,angle(Hcscheck)*1 80/pi) 
hold  on 

plot(Freqs,angle(Hstar_cssynth)*1 80/pi) 
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title(‘ Phase  \Theta_{is  is}  (\Omega)’) 
figure 

subplot(211) 

suptitle(‘Modified  ISET  FRF  response  check') 
plot(Freqs,20*log10(abs(Hischeck)),’-') 
hold  on 

plot(Freqs,20*log10(abs(Hstar_issynth)),’-.') 
legend(‘Check’, ’Synthesized’, ’Location’, ’Best’) 
if  ABC.Spring~=0 

title([‘|H_{‘,num2str(iset(end-1)),’  ‘,num2str(iset(end-1 )),... 

‘}(\Omega)|for  k=‘,num2str(ABC.Spring(iCntlter)), ’in/lb  spring’]) 

else 

title([‘|H_{‘,num2str(iset(end-1)),’  ‘,num2str(iset(end-1 )),... 
‘}(\Omega)|for  m=‘,num2str(ABC.Mass(iCntlter)),’lbm  mass’]) 
end 

xlabel('Frequency  (Hz)’) 
subplot(212) 

plot(Freqs,angle(Hischeck)*1 80/pi) 
hold  on 

plot(Freqs,angle(Hstar_issynth)*1 80/pi) 

title(‘ Phase  \Theta_{is  is}  (\Omega)’) 

fig_num=fig_num+1 ; 
end 

vars(iCntlter+2)={‘Pin’}; 

vars(iCntlter+3)={'lnital’}; 

figure(4) 

%  subplot(21 1) 

suptitle('FRF  response  change') 

%  plot(Freqs,20*log1 0(abs(HccStar_pin))) 

%  hold  on 

plot(Freqs,20*log10(abs(squeeze(Beam.H(cset,cset, :))))) 
hold  on 


o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

%  Testing  only 

axis([0  150  -160  40]) 
grid  on 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

title([‘|H_{‘,num2str(cset),’  ‘,num2str(cset),’}(\Omega)|  ]) 
xlabel('Frequency  (Hz)’) 
legend(vars(2:end), ’Location’, ’EastOutside’) 
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APPENDIX  H.  SENSITIVITY  CALCULATION  FOR  SYNTHESIZED 

BEAM 


function  [S, dLam,Pivot]=ABC_Sensistivity(Beam, ABC, LFreq, Inst, Damage, BC, Pivot) 


o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


% 


%  Variable  Description: 

%  INPUT  parameters: 

%  Beam:  Properties  of  the  beam  with  the  following  entries: 

%  phiFE  or  LimitedPhi:  Eigenvectors  of  the  Stiffness  and  Mass  matricies 
%  KGIobal:  Global  Stiffness  Matrix 

%  El:  Modulus  of  Elasticity  multiplied  by  the  crossectional 
%  moment  of  inertia 

%  length  :  length  of  the  beam 
%  width  :  Width  of  the  beam 
%  depth  :  Depth  of  the  beam 
%  elem:  Number  of  elements  in  the  beam 
%  nodes:  Array  of  nodal  locations 

%  ABC:  Properties  of  the  Artifical  Boundary  Conditions 
%  SpringMag:  The  magnitude  of  artifical  spring-to-ground  stiffness 
%  MassMag:  The  magnitude  of  artifical  dimensionless  mass 
%  Node:  Node  to  apply  the  Artifical  beam  modifications 
%  dof:Degree  of  Freedom  associated  to  the  Node  for  Beam 
%  Modifications 

%  Interval:  How  many  steps  from  1  to  the  magnitude  of 
%  modification 

%  Spring:  An  array  of  spring  values  going  from  1  to  the  magnitude 
%  of  modification  with  ABC. Interval  steps 

%  Mass:  An  array  of  mass  values  going  from  1  to  the  magnitude 
%  of  modification  with  ABC. Interval  steps 

%  Inst:  Properties  of  the  instruments  with  te  following  entries: 

%  assy_mass:  Any  additional  mass  located  at  the  nodes  of  the  beam 
%  (e.g.  Acceleromerters) 

%  massl:  Innertal  mass  of  the  additional  mass,  againlocated  at 
%  the  nodes 

%  Damage:  Properties  of  the  damage  to  the  beam 
%  elem:  Element  that  damage  occurs  on 
%  Mag:  Manitude  of  damage 
%  BC:  Boundary  Condition  applied  to  the  beam 

% 


%  OUTPUT  parameters: 

%  S:  Eigen  Sensitivity  matrixof  the  beam  before  damage.  This  is 
%  a  composite  of  every  interval  of  the  Spring  or  Mass  array 

%  dLam:  Change  in  resonance  frequencies  of  the  beam  before  and 
%  after  damage.  This  is  a  composite  of  every  interval  of  the 

%  Spring  or  Mass  array 

%  Pivot:  QR  pivot  properties  of  the  modified  beam 


o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 
/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


MGIobal^Beam.MGIobal; 
KGIobal=Beam. KGIobal; 
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S=Beam.S; 

Bmodes=ABC.Bmodes;  %  Baseline  higher  freq  modes  (+1)  to  keep 
Mmodes=ABC.Mmodes;  %  Modified  higher  freq  modes  (+1 )  to  keep 

fn_diff=zeros(length(Beam.LimitedFn),  ABC.  Interval); 
fn_rcrd=zeros((length(Beam.LimitedFn)+1),  ABC.  Interval); 

[phiTest,  lamTest,  fnTest]=damage_beam(Beam,  Inst,  BC,  Damage); 
[L_phitest,L_lamtest,L_fntest]=sys_limit(phiTest,lamTest,fnTest,LFreq); 
dLam=L_lamtest-Beam. Limited  Lam; 

Limitedphi=Beam.LimitedPhi; 

Limitediam=Beam.  Limited  Lam; 

Limitedfn=Beam.LimitedFn; 


o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

%Modify  the  Baseline  S  matrix  to  account  only  for  some  modes. 
%  %Lower  Modes 

%  S=S(2:Bmodes,:); 

%  dLam=dLam(2:Bmodes); 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

%  %  Higher  Modes 

% 

S=S(end-Bmodes:end,:); 

dLam=dLam(end-Bmodes:end); 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


for  cnt=1  :ABC. Interval 
%%  Damage  Beam 

D_mag=(Damage.mag)/1 00; 

EI_D=Beam.EI; 

EI_D(Damage.elem)=EI_D(Damage.elem)-(EI_D(Damage.elem)*D_mag); 

[Ktest,Mtest]=KandMbeam_Global(Beam,EI_D,lnst,BC); 

if  ABC.MassMag>0 
mdelta=ABC.Mass(cnt); 

Mtest(ABC.dof,ABC.dof)=Mtest(ABC.dof,ABC.dof)+mdelta; 

end 

Beam.KGIobal=KGIobal; 
if  ABC.SpringMag>0 
kdelta=ABC.Spring(cnt); 

Ktest(ABC.dof,ABC.dof)=Ktest(ABC.dof,ABC.dof)+kdelta; 

end 

%  Determine  eigenvalues  for  modified  beam,  damaged 
[phiTest_mod,lamTest_mod,fnTest_mod]  =  eig_fcn(Ktest,Mtest); 

[Lphitest_mod,Llamtest_mod,Lfntest_mod]=sys_limit(phiTest_mod,lamTest_mod,fnTest_mod,LFr 

eq); 
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%%  ABC  Conditions 
Beam.MGIobal=MGIobal; 
if  ABC.MassMag>0 
mdelta=ABC.Mass(cnt); 

Beam.  MGIobal(ABC.  dot,  ABC.  dof)=Beam.MGIobal(ABC.dof,  ABC.  dof)+mdelta; 
end 

Beam.KGIobal=KGIobal; 
if  ABC.SpringMag>0 
kdelta=ABC.Spring(cnt) 

Beam.  KGIobal(ABC.dof,  ABC. dof)=Beam.KGIobal(ABC.dof,  ABC. dof)+kdelta; 
end 

Beam.LimitedPhi=Limitedphi; 

Beam. Limited  Lam=Limitedlam; 

Beam.LimitedFn=Limitedfn; 

%  Determine  eigenvalues  for  modified  beam,  undamaged 
[phi_mod,lam_mod,fn_rnod]=eig_fcn(Beam.KGIobal,Beam.MGIobal); 

[Limitedphi_mod,Limitedlam_mod,Limitedfn_mod]=sys_limit(phi_mod,lam_mod,fn_mod,LFreq); 

Beam.LimitedPhi=Limitedphi_mod; 

Beam.  Limited  Lam=Limitedlam_mod; 

Beam.LimitedFn=Limitedfn_mod; 

dlam_Mod=Llamtest_mod(1  :size(Limitedlam_mod))-Limitedlam_mod; 
dfn(1  :length(Limitedfn_mod),cnt)=real(sqrt(dlam_Mod))/2*pi; 


%%  Sensitivity  Modification 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

fn_rcrd(:,cnt)=fn_mod(2:(length(Limitedfn)+2)); 

fn_diff(1  :length(Limitedfn_mod),cnt)=Limitedfn_mod-Lfntest_mod(1  :length(Limitedfn_mod)); 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

[S_Mod]=Sensitivity_fcn(Beam,  Inst, BC,  ABC, cnt); 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

%Modify  the  S  matrix  to  account  only  for  some  modes. 

%  %Lower  Modes 

%  S_Mod=S_Mod(2:Mmodes,:) 

%  dlam_Mod=dlam_Mod(2:Mmodes); 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

%  %  Higher  Modes 

% 

S_Mod=S_Mod(end-Mmodes:end,:); 

dlam_Mod=dlam_Mod(end-Mmodes:end); 

o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/  o/ 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


Pivot.  spring_l(cnt,:)]=qr(S_Mod,0); 
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dLam=cat(1  ,dLam,dlam_Mod); 
S=cat(1,S,S_Mod); 

Srank=rank(S); 

[~,~,Pivot.spring(cnt,:)]=qr(S,0); 


end 

display(‘Changes  in  frequency  before  and  after  damage  with  each  spring  progression  are 

ABC. Spring 

fndiff 

fn_rcrd 

end 
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